1.0 – Simulation 1: G and P Comparisons Within and Across Environments

1.1 – Setting up and Loading the Packages and Data

1.1.1 – Packages and Functions

Lets first load the relevant packages and functions. Note that, because the .Rmd file is located within the R/ folder we need to modify the source script for the “homemade” functions. Note that we have # out kableExtra as this only needs to be done once.

    rm(list=ls())
    source("./R/GP_compare_func.R")
    source("./R/evolvability_func.R")
    pacman::p_load(grid, png, shape, car, tidyr, dplyr, kableExtra, latex2exp, metafor, brms, rstan, timetree, MCMCglmm, data.table, ggplot2, knitr, rmarkdown)
    #devtools::install_github("haozhu233/kableExtra")

1.1.2 – Loading the Matrices from the data/ Folder

We want to be able to get into all the paper folders, bring in all the matrices and information about traits, perform operations on the matrices and extract the relevant information from them so we can do Monte Carlo simulations. To start, we can acquire the path names of all the matrix and description files within each paper folder of the matrices in data_GP_compare/.The below code will use a function that is custom built to grab all the above information.

      matrices <- files(dir = "./data_GP_compare/", type = "matrix", file_pattern = "[WRs][BR]")
    descripDat <- files(dir = "./data_GP_compare/", type = "D", file_pattern = "[WRs][BR]")

1.1.3 – Load the Meta-Data

The meta data is needed as it contains the relevant moderators and information for each of the papers that are included in order to merge this data with the simulations before running the meta-regression models. We can also see the data in the table below:

    meta <- read.csv("./data_GP_compare/stdy_metadata_final.csv", stringsAsFactor = FALSE)
    kable(meta) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "200px")
file authors year study env type n avg_n_env avg_n_type trait_num col genus species design_raw design group name.type nam.env env.comp
./data_GP_compare/RR073/B_n=68_G.csv Evans et al. 2015 RR073 B G 68 114.5 66.0 8 2 Poecilia reticulata half-sib half-sib fish RR073.G RR073.B BS
./data_GP_compare/RR073/B_n=161_P.csv Evans et al. 2015 RR073 B P 161 114.5 153.5 8 1 Poecilia reticulata half-sib half-sib fish RR073.P RR073.B BS
./data_GP_compare/RR073/S1_n=64_G.csv Evans et al. 2015 RR073 S G 64 105.0 66.0 8 4 Poecilia reticulata half-sib half-sib fish RR073.G RR073.S BS
./data_GP_compare/RR073/S1_n=146_P.csv Evans et al. 2015 RR073 S P 146 105.0 153.5 8 3 Poecilia reticulata half-sib half-sib fish RR073.P RR073.S BS
./data_GP_compare/RR076/B_n=69_G.csv Guan et al. 2016 RR076 B G 69 1097.0 69.0 3 6 Scophthalmus maximus half-sib half-sib fish RR076.G RR076.B BS
./data_GP_compare/RR076/B_n=2125_P.csv Guan et al. 2016 RR076 B P 2125 1097.0 2525.0 3 5 Scophthalmus maximus half-sib half-sib fish RR076.P RR076.B BS
./data_GP_compare/RR076/S1_n=69_G.csv Guan et al. 2016 RR076 S G 69 1497.0 69.0 3 8 Scophthalmus maximus half-sib half-sib fish RR076.G RR076.S BS
./data_GP_compare/RR076/S1_n=2925_P.csv Guan et al. 2016 RR076 S P 2925 1497.0 2525.0 3 7 Scophthalmus maximus half-sib half-sib fish RR076.P RR076.S BS
./data_GP_compare/RR087a/A_n=100_G.csv Sae-Lim et al. 2013 RR087a A G 100 1232.0 100.0 3 9 Oncorhynchus mykiss full-sib full-sib fish RR087a.G RR087a.A AN
./data_GP_compare/RR087a/A_n=2364_P.csv Sae-Lim et al. 2013 RR087a A P 2364 1232.0 2127.5 3 10 Oncorhynchus mykiss full-sib full-sib fish RR087a.P RR087a.A AN
./data_GP_compare/RR087a/N1_n=100_G.csv Sae-Lim et al. 2013 RR087a N G 100 995.5 100.0 3 11 Oncorhynchus mykiss full-sib full-sib fish RR087a.G RR087a.N AN
./data_GP_compare/RR087a/N1_n=1891_P.csv Sae-Lim et al. 2013 RR087a N P 1891 995.5 2127.5 3 12 Oncorhynchus mykiss full-sib full-sib fish RR087a.P RR087a.N AN
./data_GP_compare/RR087b/A_n=100_G.csv Sae-Lim et al. 2013 RR087b A G 100 1232.0 100.0 3 13 Oncorhynchus mykiss full-sib full-sib fish RR087b.G RR087b.A AN
./data_GP_compare/RR087b/A_n=2364_P.csv Sae-Lim et al. 2013 RR087b A P 2364 1232.0 2577.0 3 14 Oncorhynchus mykiss full-sib full-sib fish RR087b.P RR087b.A AN
./data_GP_compare/RR087b/N2_n=100_G.csv Sae-Lim et al. 2013 RR087b N G 100 1445.0 100.0 3 15 Oncorhynchus mykiss full-sib full-sib fish RR087b.G RR087b.N AN
./data_GP_compare/RR087b/N2_n=2790_P.csv Sae-Lim et al. 2013 RR087b N P 2790 1445.0 2577.0 3 16 Oncorhynchus mykiss full-sib full-sib fish RR087b.P RR087b.N AN
./data_GP_compare/RR087c/A_n=100_G.csv Sae-Lim et al. 2013 RR087c A G 100 1232.0 100.0 3 17 Oncorhynchus mykiss full-sib full-sib fish RR087c.G RR087c.A AN
./data_GP_compare/RR087c/A_n=2364_P.csv Sae-Lim et al. 2013 RR087c A P 2364 1232.0 2091.0 3 18 Oncorhynchus mykiss full-sib full-sib fish RR087c.P RR087c.A AN
./data_GP_compare/RR087c/N3_n=100_G.csv Sae-Lim et al. 2013 RR087c N G 100 959.0 100.0 3 19 Oncorhynchus mykiss full-sib full-sib fish RR087c.G RR087c.N AN
./data_GP_compare/RR087c/N3_n=1818_P.csv Sae-Lim et al. 2013 RR087c N P 1818 959.0 2091.0 3 20 Oncorhynchus mykiss full-sib full-sib fish RR087c.P RR087c.N AN
./data_GP_compare/sR097/B_n=26_G.csv Kasule 1991 sR097 B G 26 169.0 26.0 2 21 Dysdercus fasciatus half-sib half-sib insect sR097.G sR097.B BS
./data_GP_compare/sR097/B_n=312_P.csv Kasule 1991 sR097 B P 312 169.0 312.0 2 22 Dysdercus fasciatus half-sib half-sib insect sR097.P sR097.B BS
./data_GP_compare/sR097/S1_n=26_G.csv Kasule 1991 sR097 S G 26 169.0 26.0 2 23 Dysdercus fasciatus half-sib half-sib insect sR097.G sR097.S BS
./data_GP_compare/sR097/S1_n=312_P.csv Kasule 1991 sR097 S P 312 169.0 312.0 2 24 Dysdercus fasciatus half-sib half-sib insect sR097.P sR097.S BS
./data_GP_compare/sR110b/B_n=37_G.csv Blanckenhorn & Heyland 2004 sR110b B G 37 82.5 32.0 2 26 Scathophaga stercoraria half-sib half-sib insect sR110b.G sR110b.B BS
./data_GP_compare/sR110b/B_n=128_P.csv Blanckenhorn & Heyland 2004 sR110b B P 128 82.5 113.0 2 25 Scathophaga stercoraria half-sib half-sib insect sR110b.P sR110b.B BS
./data_GP_compare/sR110b/S1_n=27_G.csv Blanckenhorn & Heyland 2004 sR110b S G 27 62.5 32.0 2 27 Scathophaga stercoraria half-sib half-sib insect sR110b.G sR110b.S BS
./data_GP_compare/sR110b/S1_n=98_P.csv Blanckenhorn & Heyland 2004 sR110b S P 98 62.5 113.0 2 28 Scathophaga stercoraria half-sib half-sib insect sR110b.P sR110b.S BS
./data_GP_compare/sR110c/B_n=40_G.csv Blanckenhorn & Heyland 2004 sR110c B G 40 105.5 36.5 2 30 Scathophaga stercoraria half-sib half-sib insect sR110c.G sR110c.B BS
./data_GP_compare/sR110c/B_n=171_P.csv Blanckenhorn & Heyland 2004 sR110c B P 171 105.5 150.0 2 29 Scathophaga stercoraria half-sib half-sib insect sR110c.P sR110c.B BS
./data_GP_compare/sR110c/S1_n=33_G.csv Blanckenhorn & Heyland 2004 sR110c S G 33 81.0 36.5 2 32 Scathophaga stercoraria half-sib half-sib insect sR110c.G sR110c.S BS
./data_GP_compare/sR110c/S1_n=129_P.csv Blanckenhorn & Heyland 2004 sR110c S P 129 81.0 150.0 2 31 Scathophaga stercoraria half-sib half-sib insect sR110c.P sR110c.S BS
./data_GP_compare/sR115/A_n=350_G.csv Conner et al. 2003 sR115 A G 350 741.5 471.5 6 34 Raphanus raphanistrum half-sib half-sib plant sR115.G sR115.A AN
./data_GP_compare/sR115/A_n=1133_P.csv Conner et al. 2003 sR115 A P 1133 741.5 1010.5 6 33 Raphanus raphanistrum half-sib half-sib plant sR115.P sR115.A AN
./data_GP_compare/sR115/N_n=593_G.csv Conner et al. 2003 sR115 N G 593 740.5 471.5 6 35 Raphanus raphanistrum half-sib half-sib plant sR115.G sR115.N AN
./data_GP_compare/sR115/N_n=888_P.csv Conner et al. 2003 sR115 N P 888 740.5 1010.5 6 36 Raphanus raphanistrum half-sib half-sib plant sR115.P sR115.N AN
./data_GP_compare/sR116/A_n=366_G.csv Delcourt et al. 2010 sR116 A G 366 1840.5 365.5 8 38 Drosophila serrata half-sib half-sib insect sR116.G sR116.A AN
./data_GP_compare/sR116/A_n=3315_P.csv Delcourt et al. 2010 sR116 A P 3315 1840.5 3299.0 8 37 Drosophila serrata half-sib half-sib insect sR116.P sR116.A AN
./data_GP_compare/sR116/N1_n=365_G.csv Delcourt et al. 2010 sR116 N G 365 1824.0 365.5 8 40 Drosophila serrata half-sib half-sib insect sR116.G sR116.N AN
./data_GP_compare/sR116/N1_n=3283_P.csv Delcourt et al. 2010 sR116 N P 3283 1824.0 3299.0 8 39 Drosophila serrata half-sib half-sib insect sR116.P sR116.N AN
./data_GP_compare/sR145/B_n=52_G.csv Tucic et al. 1991 sR145 B G 52 101.0 47.5 8 42 Acanthoscelides obtectus full-sib full-sib insect sR145.G sR145.B BS
./data_GP_compare/sR145/B_n=150_P.csv Tucic et al. 1991 sR145 B P 150 101.0 137.0 8 41 Acanthoscelides obtectus full-sib full-sib insect sR145.P sR145.B BS
./data_GP_compare/sR145/S1_n=43_G.csv Tucic et al. 1991 sR145 S G 43 83.5 47.5 8 44 Acanthoscelides obtectus full-sib full-sib insect sR145.G sR145.S BS
./data_GP_compare/sR145/S1_n=124_P.csv Tucic et al. 1991 sR145 S P 124 83.5 137.0 8 43 Acanthoscelides obtectus full-sib full-sib insect sR145.P sR145.S BS
./data_GP_compare/WB002a/A_n=39_G.csv Begin & Roff 2001 WB002a A G 39 372.0 39.0 5 45 Gryllus pennsylvanicus full-sib full-sib insect WB002a.G WB002a.A AN
./data_GP_compare/WB002a/A_n=705_P.csv Begin & Roff 2001 WB002a A P 705 372.0 605.0 5 46 Gryllus pennsylvanicus full-sib full-sib insect WB002a.P WB002a.A AN
./data_GP_compare/WB002a/N1_n=39_G.csv Begin & Roff 2001 WB002a N G 39 272.0 39.0 5 47 Gryllus pennsylvanicus full-sib full-sib insect WB002a.G WB002a.N AN
./data_GP_compare/WB002a/N1_n=505_P.csv Begin & Roff 2001 WB002a N P 505 272.0 605.0 5 48 Gryllus pennsylvanicus full-sib full-sib insect WB002a.P WB002a.N AN
./data_GP_compare/WB003/A_n=56_G.csv Begin et al. 2004 WB003 A G 56 294.5 59.0 5 50 Gryllus firmus full-sib full-sib insect WB003.G WB003.A AN
./data_GP_compare/WB003/A_n=533_P.csv Begin et al. 2004 WB003 A P 533 294.5 697.5 5 49 Gryllus firmus full-sib full-sib insect WB003.P WB003.A AN
./data_GP_compare/WB003/N1_n=62_G.csv Begin et al. 2004 WB003 N G 62 387.5 59.0 5 51 Gryllus firmus full-sib full-sib insect WB003.G WB003.N AN
./data_GP_compare/WB003/N1_n=862_P.csv Begin et al. 2004 WB003 N P 862 387.5 697.5 5 52 Gryllus firmus full-sib full-sib insect WB003.P WB003.N AN
./data_GP_compare/WB003_2/A_n=56_G.csv Begin et al. 2004 WB003_2 A G 56 294.5 58.0 5 54 Gryllus firmus full-sib full-sib insect WB003_2.G WB003_2.A AN
./data_GP_compare/WB003_2/A_n=533_P.csv Begin et al. 2004 WB003_2 A P 533 294.5 624.0 5 53 Gryllus firmus full-sib full-sib insect WB003_2.P WB003_2.A AN
./data_GP_compare/WB003_2/N2_n=60_G.csv Begin et al. 2004 WB003_2 N G 60 462.0 58.0 5 55 Gryllus firmus full-sib full-sib insect WB003_2.G WB003_2.N AN
./data_GP_compare/WB003_2/N2_n=715_P.csv Begin et al. 2004 WB003_2 N P 715 462.0 624.0 5 56 Gryllus firmus full-sib full-sib insect WB003_2.P WB003_2.N AN
./data_GP_compare/WB003a/A_n=45_G.csv Begin et al. 2004 WB003a A G 45 153.5 43.0 5 58 Gryllus firmus full-sib full-sib insect WB003a.G WB003a.A AN
./data_GP_compare/WB003a/A_n=262_P.csv Begin et al. 2004 WB003a A P 262 153.5 288.5 5 57 Gryllus firmus full-sib full-sib insect WB003a.P WB003a.A AN
./data_GP_compare/WB003a/N1_n=41_G.csv Begin et al. 2004 WB003a N G 41 82.5 43.0 5 60 Gryllus firmus full-sib full-sib insect WB003a.G WB003a.N AN
./data_GP_compare/WB003a/N1_n=315_P.csv Begin et al. 2004 WB003a N P 315 82.5 288.5 5 59 Gryllus firmus full-sib full-sib insect WB003a.P WB003a.N AN
./data_GP_compare/WB003a_2/A_n=45_G.csv Begin et al. 2004 WB003a_2 A G 45 153.5 38.5 5 62 Gryllus firmus full-sib full-sib insect WB003a_2.G WB003a_2.A AN
./data_GP_compare/WB003a_2/A_n=262_P.csv Begin et al. 2004 WB003a_2 A P 262 153.5 197.5 5 61 Gryllus firmus full-sib full-sib insect WB003a_2.P WB003a_2.A AN
./data_GP_compare/WB003a_2/N2_n=32_G.csv Begin et al. 2004 WB003a_2 N G 32 178.0 38.5 5 64 Gryllus firmus full-sib full-sib insect WB003a_2.G WB003a_2.N AN
./data_GP_compare/WB003a_2/N2_n=133_P.csv Begin et al. 2004 WB003a_2 N P 133 178.0 197.5 5 63 Gryllus firmus full-sib full-sib insect WB003a_2.P WB003a_2.N AN
./data_GP_compare/WB004/A_n=300_G.csv Bennington & McGraw 2000 WB004 A G 300 342.0 300.0 5 65 Impatiens pallida half-sib half-sib plant WB004.G WB004.A AN
./data_GP_compare/WB004/A_n=384_P.csv Bennington & McGraw 2000 WB004 A P 384 342.0 244.5 5 66 Impatiens pallida half-sib half-sib plant WB004.P WB004.A AN
./data_GP_compare/WB004/N1_n=300_G.csv Bennington & McGraw 2000 WB004 N G 300 202.5 300.0 5 68 Impatiens pallida half-sib half-sib plant WB004.G WB004.N AN
./data_GP_compare/WB004/N1_n=105_P.csv Bennington & McGraw 2000 WB004 N P 105 202.5 244.5 5 67 Impatiens pallida half-sib half-sib plant WB004.P WB004.N AN
./data_GP_compare/WB004a/A_n=300_G.csv Bennington & McGraw 2000 WB004a A G 300 190.0 300.0 5 69 Impatiens pallida half-sib half-sib plant WB004a.G WB004a.A AN
./data_GP_compare/WB004a/A_n=80_P.csv Bennington & McGraw 2000 WB004a A P 80 190.0 108.5 5 70 Impatiens pallida half-sib half-sib plant WB004a.P WB004a.A AN
./data_GP_compare/WB004a/N1_n=300_G.csv Bennington & McGraw 2000 WB004a N G 300 218.5 300.0 5 72 Impatiens pallida half-sib half-sib plant WB004a.G WB004a.N AN
./data_GP_compare/WB004a/N1_n=137_P.csv Bennington & McGraw 2000 WB004a N P 137 218.5 108.5 5 71 Impatiens pallida half-sib half-sib plant WB004a.P WB004a.N AN
./data_GP_compare/WB007/B_n=132_G.csv Brock et al. 2010 WB007 B G 132 366.5 131.5 11 73 Brassica rapa RIL half-sib plant WB007.G WB007.B BS
./data_GP_compare/WB007/B_n=601_P.csv Brock et al. 2010 WB007 B P 601 366.5 562.5 11 74 Brassica rapa RIL half-sib plant WB007.P WB007.B BS
./data_GP_compare/WB007/S1_n=131_G.csv Brock et al. 2010 WB007 S G 131 327.5 131.5 11 75 Brassica rapa RIL half-sib plant WB007.G WB007.S BS
./data_GP_compare/WB007/S1_n=524_P.csv Brock et al. 2010 WB007 S P 524 327.5 562.5 11 76 Brassica rapa RIL half-sib plant WB007.P WB007.S BS
./data_GP_compare/WB010/A_n=141_G.csv Collins et al. 1999 WB010 A G 141 840.5 101.5 4 77 Achroia grisella half-sib half-sib insect WB010.G WB010.A AN
./data_GP_compare/WB010/A_n=1540_P.csv Collins et al. 1999 WB010 A P 1540 840.5 991.5 4 78 Achroia grisella half-sib half-sib insect WB010.P WB010.A AN
./data_GP_compare/WB010/N1_n=62_G.csv Collins et al. 1999 WB010 N G 62 252.5 101.5 4 80 Achroia grisella half-sib half-sib insect WB010.G WB010.N AN
./data_GP_compare/WB010/N1_n=443_P.csv Collins et al. 1999 WB010 N P 443 252.5 991.5 4 79 Achroia grisella half-sib half-sib insect WB010.P WB010.N AN
./data_GP_compare/WB011/B_n=404_G.csv Czesak & Fox 2003 WB011 B G 404 1535.5 404.0 2 82 Stator limbatus half-sib half-sib insect WB011.G WB011.B BS
./data_GP_compare/WB011/B_n=2667_P.csv Czesak & Fox 2003 WB011 B P 2667 1535.5 2667.0 2 81 Stator limbatus half-sib half-sib insect WB011.P WB011.B BS
./data_GP_compare/WB011/S1_n=404_G.csv Czesak & Fox 2003 WB011 S G 404 1535.5 404.0 2 84 Stator limbatus half-sib half-sib insect WB011.G WB011.S BS
./data_GP_compare/WB011/S1_n=2667_P.csv Czesak & Fox 2003 WB011 S P 2667 1535.5 2667.0 2 83 Stator limbatus half-sib half-sib insect WB011.P WB011.S BS
./data_GP_compare/WB028/B_n=60_G.csv Klause & Morin 2001 WB028 B G 60 330.0 60.0 2 85 Priophorus pallipes full-sib full-sib insect WB028.G WB028.B BS
./data_GP_compare/WB028/B_n=600_P.csv Klause & Morin 2001 WB028 B P 600 330.0 600.0 2 86 Priophorus pallipes full-sib full-sib insect WB028.P WB028.B BS
./data_GP_compare/WB028/S1_n=60_G.csv Klause & Morin 2001 WB028 S G 60 330.0 60.0 2 87 Priophorus pallipes full-sib full-sib insect WB028.G WB028.S BS
./data_GP_compare/WB028/S1_n=600_P.csv Klause & Morin 2001 WB028 S P 600 330.0 600.0 2 88 Priophorus pallipes full-sib full-sib insect WB028.P WB028.S BS
./data_GP_compare/WB029/B_n=63_G.csv King et al. 2011 WB029 B G 63 862.5 63.0 3 90 Gryllus firmus half-sib half-sib insect WB029.G WB029.B BS
./data_GP_compare/WB029/B_n=1662_P.csv King et al. 2011 WB029 B P 1662 862.5 1695.0 3 89 Gryllus firmus half-sib half-sib insect WB029.P WB029.B BS
./data_GP_compare/WB029/S1_n=63_G.csv King et al. 2011 WB029 S G 63 895.5 63.0 3 92 Gryllus firmus half-sib half-sib insect WB029.G WB029.S BS
./data_GP_compare/WB029/S1_n=1728_P.csv King et al. 2011 WB029 S P 1728 895.5 1695.0 3 91 Gryllus firmus half-sib half-sib insect WB029.P WB029.S BS
./data_GP_compare/WB033/A_n=60_G.csv Lau et al. 2014 WB033 A G 60 293.0 60.0 4 94 Arabidopsis thaliana RIL half-sib plant WB033.G WB033.A AN
./data_GP_compare/WB033/A_n=526_P.csv Lau et al. 2014 WB033 A P 526 293.0 513.0 4 93 Arabidopsis thaliana RIL half-sib plant WB033.P WB033.A AN
./data_GP_compare/WB033/N1_n=60_G.csv Lau et al. 2014 WB033 N G 60 280.0 60.0 4 96 Arabidopsis thaliana RIL half-sib plant WB033.G WB033.N AN
./data_GP_compare/WB033/N1_n=500_P.csv Lau et al. 2014 WB033 N P 500 280.0 513.0 4 95 Arabidopsis thaliana RIL half-sib plant WB033.P WB033.N AN
./data_GP_compare/WB033a/A_n=60_G.csv Lau et al. 2014 WB033a A G 60 259.5 60.0 4 98 Arabidopsis thaliana RIL half-sib plant WB033a.G WB033a.A AN
./data_GP_compare/WB033a/A_n=459_P.csv Lau et al. 2014 WB033a A P 459 259.5 464.5 4 97 Arabidopsis thaliana RIL half-sib plant WB033a.P WB033a.A AN
./data_GP_compare/WB033a/N1_n=60_G.csv Lau et al. 2014 WB033a N G 60 265.0 60.0 4 100 Arabidopsis thaliana RIL half-sib plant WB033a.G WB033a.N AN
./data_GP_compare/WB033a/N1_n=470_P.csv Lau et al. 2014 WB033a N P 470 265.0 464.5 4 99 Arabidopsis thaliana RIL half-sib plant WB033a.P WB033a.N AN
./data_GP_compare/WB035/B_n=432_G.csv Messina et al. 2013 WB035 B G 432 1316.5 432.0 3 102 Callosobruchus maculatus half-sib half-sib insect WB035.G WB035.B BS
./data_GP_compare/WB035/B_n=2201_P.csv Messina et al. 2013 WB035 B P 2201 1316.5 2204.0 3 101 Callosobruchus maculatus half-sib half-sib insect WB035.P WB035.B BS
./data_GP_compare/WB035/S1_n=432_G.csv Messina et al. 2013 WB035 S G 432 1319.5 432.0 3 104 Callosobruchus maculatus half-sib half-sib insect WB035.G WB035.S BS
./data_GP_compare/WB035/S1_n=2207_P.csv Messina et al. 2013 WB035 S P 2207 1319.5 2204.0 3 103 Callosobruchus maculatus half-sib half-sib insect WB035.P WB035.S BS
./data_GP_compare/WB038/B_n=135_G.csv Rauter&Moore 2002 WB038 B G 135 387.0 135.0 4 105 Nicrophorus pustulatus half-sib half-sib insect WB038.G WB038.B BS
./data_GP_compare/WB038/B_n=639_P.csv Rauter&Moore 2002 WB038 B P 639 387.0 639.0 4 106 Nicrophorus pustulatus half-sib half-sib insect WB038.P WB038.B BS
./data_GP_compare/WB038/S1_n=135_G.csv Rauter&Moore 2002 WB038 S G 135 387.0 135.0 4 107 Nicrophorus pustulatus half-sib half-sib insect WB038.G WB038.S BS
./data_GP_compare/WB038/S1_n=636_P.csv Rauter&Moore 2002 WB038 S P 636 387.0 639.0 4 108 Nicrophorus pustulatus half-sib half-sib insect WB038.P WB038.S BS
./data_GP_compare/WB045/A_n=39_G.csv Simons & Roff 1996 WB045 A G 39 209.5 39.0 7 110 Gryllus pennsylvanicus full-sib full-sib insect WB045.G WB045.A AN
./data_GP_compare/WB045/A_n=380_P.csv Simons & Roff 1996 WB045 A P 380 209.5 371.0 7 109 Gryllus pennsylvanicus full-sib full-sib insect WB045.P WB045.A AN
./data_GP_compare/WB045/N1_n=39_G.csv Simons & Roff 1996 WB045 N G 39 200.5 39.0 7 112 Gryllus pennsylvanicus full-sib full-sib insect WB045.G WB045.N AN
./data_GP_compare/WB045/N1_n=362_P.csv Simons & Roff 1996 WB045 N P 362 200.5 371.0 7 111 Gryllus pennsylvanicus full-sib full-sib insect WB045.P WB045.N AN
./data_GP_compare/WB045a/A_n=39_G.csv Simons & Roff 1996 WB045a A G 39 215.0 39.0 6 113 Gryllus pennsylvanicus full-sib full-sib insect WB045a.G WB045a.A AN
./data_GP_compare/WB045a/A_n=391_P.csv Simons & Roff 1996 WB045a A P 391 215.0 381.0 6 114 Gryllus pennsylvanicus full-sib full-sib insect WB045a.P WB045a.A AN
./data_GP_compare/WB045a/N1_n=39_G.csv Simons & Roff 1996 WB045a N G 39 205.0 39.0 6 116 Gryllus pennsylvanicus full-sib full-sib insect WB045a.G WB045a.N AN
./data_GP_compare/WB045a/N1_n=371_P.csv Simons & Roff 1996 WB045a N P 371 205.0 381.0 6 115 Gryllus pennsylvanicus full-sib full-sib insect WB045a.P WB045a.N AN
./data_GP_compare/WB053/A_n=24_G.csv Via & Conner 1995 WB053 A G 24 252.0 24.0 2 117 Tribolium castaneum half-sib half-sib insect WB053.G WB053.A AN
./data_GP_compare/WB053/A_n=480_P.csv Via & Conner 1995 WB053 A P 480 252.0 480.0 2 118 Tribolium castaneum half-sib half-sib insect WB053.P WB053.A AN
./data_GP_compare/WB053/N1_n=24_G.csv Via & Conner 1995 WB053 N G 24 252.0 24.0 2 119 Tribolium castaneum half-sib half-sib insect WB053.G WB053.N AN
./data_GP_compare/WB053/N1_n=480_P.csv Via & Conner 1995 WB053 N P 480 252.0 480.0 2 120 Tribolium castaneum half-sib half-sib insect WB053.P WB053.N AN
./data_GP_compare/WB053_3/A_n=24_G.csv Via & Conner 1995 WB053_3 A G 24 252.0 24.0 2 121 Tribolium castaneum half-sib half-sib insect WB053_3.G WB053_3.A AN
./data_GP_compare/WB053_3/A_n=480_P.csv Via & Conner 1995 WB053_3 A P 480 252.0 480.0 2 122 Tribolium castaneum half-sib half-sib insect WB053_3.P WB053_3.A AN
./data_GP_compare/WB053_3/N3_n=24_G.csv Via & Conner 1995 WB053_3 N G 24 252.0 24.0 2 123 Tribolium castaneum half-sib half-sib insect WB053_3.G WB053_3.N AN
./data_GP_compare/WB053_3/N3_n=480_P.csv Via & Conner 1995 WB053_3 N P 480 252.0 480.0 2 124 Tribolium castaneum half-sib half-sib insect WB053_3.P WB053_3.N AN
./data_GP_compare/WB053b/A_n=80_G.csv Via & Conner 1995 WB053b A G 80 280.0 80.0 2 126 Tribolium castaneum half-sib half-sib insect WB053b.G WB053b.A AN
./data_GP_compare/WB053b/A_n=480_P.csv Via & Conner 1995 WB053b A P 480 280.0 480.0 2 125 Tribolium castaneum half-sib half-sib insect WB053b.P WB053b.A AN
./data_GP_compare/WB053b/N1_n=80_G.csv Via & Conner 1995 WB053b N G 80 280.0 80.0 2 128 Tribolium castaneum half-sib half-sib insect WB053b.G WB053b.N AN
./data_GP_compare/WB053b/N1_n=480_P.csv Via & Conner 1995 WB053b N P 480 280.0 480.0 2 127 Tribolium castaneum half-sib half-sib insect WB053b.P WB053b.N AN
./data_GP_compare/WB053b_3/A_n=80_G.csv Via & Conner 1995 WB053b_3 A G 80 280.0 80.0 2 130 Tribolium castaneum half-sib half-sib insect WB053b_3.G WB053b_3.A AN
./data_GP_compare/WB053b_3/A_n=480_P.csv Via & Conner 1995 WB053b_3 A P 480 280.0 480.0 2 129 Tribolium castaneum half-sib half-sib insect WB053b_3.P WB053b_3.A AN
./data_GP_compare/WB053b_3/N3_n=80_G.csv Via & Conner 1995 WB053b_3 N G 80 280.0 80.0 2 132 Tribolium castaneum half-sib half-sib insect WB053b_3.G WB053b_3.N AN
./data_GP_compare/WB053b_3/N3_n=480_P.csv Via & Conner 1995 WB053b_3 N P 480 280.0 480.0 2 131 Tribolium castaneum half-sib half-sib insect WB053b_3.P WB053b_3.N AN
./data_GP_compare/WB057/A_n=22_G.csv Grill et al. 1997 WB057 A G 22 149.5 24.0 4 133 Harmonia axyridis full-sib full-sib insect WB057.G WB057.A AN
./data_GP_compare/WB057/A_n=277_P.csv Grill et al. 1997 WB057 A P 277 149.5 277.0 4 134 Harmonia axyridis full-sib full-sib insect WB057.P WB057.A AN
./data_GP_compare/WB057/N1_n=26_G.csv Grill et al. 1997 WB057 N G 26 151.5 24.0 4 135 Harmonia axyridis full-sib full-sib insect WB057.G WB057.N AN
./data_GP_compare/WB057/N1_n=277_P.csv Grill et al. 1997 WB057 N P 277 151.5 277.0 4 136 Harmonia axyridis full-sib full-sib insect WB057.P WB057.N AN

file – Described the path to the filename for the specific matrix that has been extracted from the paper.

authors – First name of author for the paper matrices were extracted from.

year – Year study was published

study – Study identifier

env – Environment traits were measured. B = Benign; S = stressful; A = non-novel; N = novel. Note that B and S are really just more detailed descriptors of the environments being stressful based on reading of the paper and the authors’ understanding of the system.

type – The type of matrix.G = Genetic variance/covariance matrix; P Phenotypic variance/covariance matrix

n – The average sample size. For P this is the total number of individuals, if G this is the number of families.

avg_n_env – Is the average sample size across environments

avg_n_type – Is the average sample size across the G and P matrices

trait_num – The number of traits in the matrices

col – Column indicator which matches names of matrices with meta-data.

genus Genus name

species Binomial species name

design_raw – The raw breeding design described in the paper. Can be half-sib, full-sib or RIL (Recombinant Inbred Lines). These are recoded so that RIL is just part of a half-sib design

design – Simplified breeding design

group – General taxonomic group studied

name.type – Interaction/ composite variable of study number and matrix type.

name.env – Interaction/ composite variable of study number and environment type.

env.comp – Environmental comparison moderator; BS = benign/stressful comparisons; AN = non-novel/novel comparison.

1.1.4 – Importing the Matrices

Import matrices and relevant data.

         mats <- import_mat(matrices)
        dscrp <- import_mat(descripDat)

Scale the data to prevent disproportionately large effects of the traits with large values. Not scaling will likely lead to similar results, but this is a safer option.

  standard <- list()
  for (i in 1:length(dscrp)){
    # standardize descriptive data
    # calculate the average mean for each trait
    dscrp[[i]]$stanfac <- ave(dscrp[[i]]$mean,dscrp[[i]]$trait)
    #standardize mean, sd, se, Va
    dscrp[[i]]$mean <- dscrp[[i]]$mean/dscrp[[i]]$stanfac
    dscrp[[i]]$sd <- dscrp[[i]]$sd/dscrp[[i]]$stanfac
    dscrp[[i]]$se <- dscrp[[i]]$se/dscrp[[i]]$stanfac
    dscrp[[i]]$Va <- dscrp[[i]]$Va/(dscrp[[i]]$stanfac^2)
    # save the standardization factor in to another list
    standard[[i]] <- data.frame(trait=unique(dscrp[[i]]$trait),stanfac=dscrp[[i]]$stanfac[which(!duplicated(dscrp[[i]]$trait))]) 
    # remove spaces and dashes from the trait names to be consistent with the names in the covariance-matrices.
    standard[[i]]$trait <- gsub(" ","\\.", as.character(standard[[i]]$trait))
    standard[[i]]$trait <- gsub("-","\\.", standard[[i]]$trait)
    rownames(standard[[i]]) <- standard[[i]]$trait
    }
  
  id_d <- sapply(names(dscrp),function(x){strsplit(x,"/")[[1]][3]})
  id_m <- sapply(names(mats),function(x){strsplit(x,"/")[[1]][3]})
  for (i in 1:length(mats)){
    # select the right standarization factors
    j <- which(id_d == id_m[[i]])
    std <- standard[[j]]
    # remove traits not represented in the covar-matrix
    std <- std[which(std$trait %in% names(mats[[i]])),]
    # order traits as in the covar-matrix
    std <- std[names(mats[[i]]),]
    # check whether the standardization factors line up
    if(sum(names(mats[[i]]) != std$trait)>0){stop("There is a inconsistency between trait names of a matrix and a description file")}
    mats[[i]] <- mats[[i]]/(std$stanfac %*% t(std$stanfac))
    }
  
  # remove studies without means which therefore cannot be standardized 
  rm.ids <- which(is.na(lapply(mats,sum)))
  mats <- mats[-rm.ids]
  dscrp <- dscrp[-which(id_d %in% id_m[rm.ids])]

Extract the sample size and get a sense of the number of traits from the filenames of each list.

                n <- match_n(names(mats))                       
        trait_num <- unlist(lapply(mats, function(x) nrow(x)))

1.2 – Simulations using G and P matrices with given sample sizes

Now that we have loaded all the data we can run the simulations described in the paper. This is done iteratively for each matrix and for each effect size. This generally does not take too long, even for 5000 simulations.

        sims = 5000

        system.time(
            # Proportion of variance of the dominant eigen value
        tota_var <- mapply(function(x,y) sim(matrix = as.matrix(x), n = y, type = "totalVar", cor = FALSE, sims = sims, matrix_fix = "bend"), x = mats, y = split(n, 1:length(n))) )
##    user  system elapsed 
## 407.424   9.352 420.286
        system.time(
            # Proportion of variance of the dominant eigen value
        var_along_logit <- mapply(function(x,y) sim(matrix = as.matrix(x), n = y, type = "varAlong", cor = FALSE, logit = TRUE, sims = sims, matrix_fix = "bend"), x = mats, y = split(n, 1:length(n))) )
##    user  system elapsed 
## 554.818   9.704 568.664
        system.time(
            # raw simulated matrices 
        matrices <- mapply(function(x,y) sim(matrix = as.matrix(x), n = y, type = "raw_matrix", cor = FALSE, sims = sims, matrix_fix = "bend"), x = mats, y = split(n, 1:length(n))) )
##    user  system elapsed 
## 200.958   7.285 208.637

1.3 – Generating Effect Sizes and Sampling Variances for Meta-analysis

Now that we have got all the simulations of the G and P matrices for the three different effect sizes of interest, we can extract from these distributions the average effect size and corresponding sampling error (i.e., the standard deviation). Note that, to generate some of the effect sizes we need to rely on the raw simulated matrices (i.e., the matrices object)

First, lets create the datasets for the first two effect sizes as they are already done for us. We will plot these data to to have a look and make sure nothing crazy is happening. We will also calculate the weights for the effect sizes and add these into the data frame.

First for standardised mean difference (SMD) in variance:

# Extract the relevant columns from the metadata file
     meta_sub_env <- meta[,c("authors", "year", "study",  "env", "avg_n_env", "trait_num", "genus", "species", "design", "group", "env.comp")]

    meta_sub_type <- meta[,c("authors", "year", "study", "type",  "avg_n_type", "trait_num", "genus", "species", "design", "group", "env.comp")]

We can have a look at the simulated distributions for each study graphically and print out a list used for comparisons to make sure things are working as expected.

    create_data(tota_var, 
                type = "SMD", 
                environment = "across", 
                plot = "yes")

Figure 1.3.1 – Effect size distributions of the standardised mean differences in total variance between matrices.

Now that this looks good we can create the data, add important information to the data frame, and save the datafile so that we no longer need to re-run these simulations each time we want to re-run analyses.

# Final data 
# Standardized mean difference in variance between G-G and P-P across environments
    
               SMD <- create_data(tota_var, 
                                  type = "SMD", 
                                  environment = "across", 
                                  plot = "no")
    
          SMD_data <- combine_meta_data(SMD, 
                                        meta_sub_type, 
                                        environment = "across")

    SMD_data$err   <- 1:nrow(SMD_data)
     SMD_data$wi   <- 1/SMD_data$S_var
    
    write.csv(SMD_data, file = "./data_GP_compare/SMD_data_standardised.csv")

Directionality of Effects: It’s important to remember here that positive values mean that the total variance in \(G_{novel}\) or \(P_{novel}\) is larger than \(G_{non-novel}\) or \(P_{non-novel}\), where as negative values mean that the total variance in \(G_{novel}\) or \(P_{novel}\) is smaller than \(G_{non-novel}\) or \(P_{non-novel}\).

Second, we can get the logit transformed proportion of variance along \(g_{max}\). This is to make the differences at the extremes (i.e., 0 and 1) less crazy compared to other values as it stretches out the ends of the distribution and normalises things as best as possible. Again, we can plot the simulations out to make sure that they look good before creating the data.

             create_data(var_along_logit, 
                        type = "raw_diff", 
                        environment = "across", 
                        plot = "yes")

Figure 1.3.2 – Effect size distributions of the logit transformed proportion of variance along \(g/p_{max}\) between matrices.

Once we are satisfied, we can create the data frame so we don’t need to run the simulations again.

          var_along_logit_d <- create_data(var_along_logit, 
                                           type = "raw_diff", 
                                           environment = "across", 
                                           plot = "no")
  
        var_along_logit_data <- combine_meta_data(var_along_logit_d, 
                                                   meta_sub_type, 
                                                   environment = "across")

    var_along_logit_data$err <- 1:nrow(var_along_logit_data)
     var_along_logit_data$wi <- 1/var_along_logit_data$S_var

    write.csv(var_along_logit_data, file = "./data_GP_compare/var_along_logit_data_standardised.csv")

Directionality of Effects: It’s important to remember here that positive values mean that the proportion of variance along \(gmax_{novel}\) or \(pmax_{novel}\) is larger than \(gmax_{non-novel}\) or \(pmax_{non-novel}\), where as negative values mean that the total variance in \(gmax_{novel}\) or \(pmax_{novel}\) is smaller than \(gmax_{non-novel}\) or \(pmax_{non-novel}\). We also need to remember that these are logit transformed.

Last, we can then calculate the angle between G-G across environments and G-P within each of the environments.

Plot to have a look at the effect size distribution across studies:

            create_data(matrices, 
                        type = "angle", 
                        environment = "across", 
                        plot = "yes")

Figure 1.3.3 – Effect size distributions of the logit transformed angle between \(g_{max}\) across environments.

The distributions are not as clean as the others, but not terrible and we can create the data for comparisons of the same matrix (e.g., \(G_{novel}\) compared to \(G_{non-novel}\)) across environments.

# Angle between P-P and G-G across environments for each study.
    
             angle_across <- create_data(matrices, 
                                            type = "angle", 
                                            environment = "across", 
                                            plot = "no")
    
        angle_across_data <- combine_meta_data(angle_across, 
                                                meta_sub_type, 
                                                environment = "across")

    angle_across_data$err <- 1:nrow(angle_across_data)
     angle_across_data$wi <- 1/angle_across_data$S_var

    write.csv(angle_across_data, file = "./data_GP_compare/angle_across_data_standardised.csv")

Directionality of Effects: Angles are logit transformed by log transformation of the angle itself and dividing this by 90 degrees. Positive effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(gmax_{novel}\) or \(pmax_{non-novel}\) and \(pmax_{novel}\) is approaching 90 degrees (i.e., values around 7 = 90 degrees). In contrast, negative effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(gmax_{novel}\) or \(pmax_{non-novel}\) and \(pmax_{novel}\) is approaching 0 degrees (i.e., values around -7 = 0 degrees).

We can then do the same for comparisons of different matrices within each environment to understand how G and P align, for example (i.e., \(G_{novel}\) compared to \(P_{novel}\))

        create_data(matrices, type = "angle", environment = "within", plot = "yes")

Figure 1.3.4 – Effect size distributions of the logit transformed angle between \(g_{max}\) and \(p_{max}\) within environments.

Once we are happy with this we can create the data frame so we again don’t need to re-run these simulations again.

    # Angle between P and G within an environment for each study.

             angle_within <- create_data(matrices, type = "angle", environment = "within", plot = "no")
    
        angle_within_data <- combine_meta_data(angle_within, meta_sub_env, environment = "within")

    angle_within_data$err <- 1:nrow(angle_within_data)
     angle_within_data$wi <- 1/angle_within_data$S_var

    write.csv(angle_within_data, file = "./data_GP_compare/angle_within_data_standardised.csv")

Directionality of Effects: Angles are logit transformed by log transformation of the angle itself and dividing this by 90 degrees. Positive effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(pmax_{non-novel}\) or \(gmax_{novel}\) and \(pmax_{novel}\) is approaching 90 degrees (i.e., values around 7 = 90 degrees). In contrast, negative effect sizes indicate that the angle between \(gmax_{non-novel}\) and \(pmax_{non-novel}\) or \(gmax_{novel}\) and \(pmax_{novel}\) is approaching 0 degrees (i.e., values around -7 = 0 degrees). Values in between are varying angles.

2.0 – Simulation 2: Alignment of Plastic Responses with Gnovel and Gnon-novel

Now that we have simulated all the data that is relevant for comparisons between G and P we can begin to ask questions about the evolvability of plastic responses. Since we only need the mean trait changes for traits across environments to generate a phenotype response vector we can actually increase the total sample size of studies included because we don’t need to rely on getting the entire P matrix for simulations.

2.1 – Importing the data

To import the data we have separated out all of the matrices for each of the two simulations and put these in separate folders. This just makes it a little more clear what data was used for each of the simulations described.

      matrices_evolve <- files(dir = "./data_evolvability/", type = "matrix", file_pattern = "[WRs][BR]")
    descripDat_evolve <- files(dir = "./data_evolvability/", type = "D", file_pattern = "[WRs][BR]")

We also need to bring in some new meta-data that includes the additional studies that are added. One can view the data below.

    meta_evolve <- read.csv("./data_evolvability/metadata_final.csv", stringsAsFactor = FALSE)
    kable(meta_evolve) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "200px")
file authors year study env type n trait_num col genus species design_raw design group env.comp name.type name.env
./data_evolvability/RR062/A_n=10_G.csv Beacham 1988 RR062 A G 10 3 1 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062.G RR062.A
./data_evolvability/RR062/A_n=100_P.csv Beacham 1988 RR062 A P 100 3 2 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062.P RR062.A
./data_evolvability/RR062/N1_n=10_G.csv Beacham 1988 RR062 N G 10 3 3 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062.G RR062.N
./data_evolvability/RR062/N1_n=74_P.csv Beacham 1988 RR062 N P 74 3 4 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062.P RR062.N
./data_evolvability/RR062_b/A_n=10_G.csv Beacham 1988 RR062_b A G 10 3 5 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062_b.G RR062_b.A
./data_evolvability/RR062_b/A_n=100_P.csv Beacham 1988 RR062_b A P 100 3 6 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062_b.P RR062_b.A
./data_evolvability/RR062_b/N2_n=10_G.csv Beacham 1988 RR062_b N G 10 3 7 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062_b.G RR062_b.N
./data_evolvability/RR062_b/N2_n=100_P.csv Beacham 1988 RR062_b N P 100 3 8 Oncorhynchus gorbuscha half-sib half-sib fish AN RR062_b.P RR062_b.N
./data_evolvability/RR062b/A_n=10_G.csv Beacham 1988 RR062b A G 10 3 9 Oncorhynchus keta half-sib half-sib fish AN RR062b.G RR062b.A
./data_evolvability/RR062b/A_n=100_P.csv Beacham 1988 RR062b A P 100 3 10 Oncorhynchus keta half-sib half-sib fish AN RR062b.P RR062b.A
./data_evolvability/RR062b/N1_n=10_G.csv Beacham 1988 RR062b N G 10 3 11 Oncorhynchus keta half-sib half-sib fish AN RR062b.G RR062b.N
./data_evolvability/RR062b/N1_n=98_P.csv Beacham 1988 RR062b N P 98 3 12 Oncorhynchus keta half-sib half-sib fish AN RR062b.P RR062b.N
./data_evolvability/RR062b_b/A_n=10_G.csv Beacham 1988 RR062b_b A G 10 3 13 Oncorhynchus keta half-sib half-sib fish AN RR062b_b.G RR062b_b.A
./data_evolvability/RR062b_b/A_n=100_P.csv Beacham 1988 RR062b_b A P 100 3 14 Oncorhynchus keta half-sib half-sib fish AN RR062b_b.P RR062b_b.A
./data_evolvability/RR062b_b/N2_n=10_G.csv Beacham 1988 RR062b_b N G 10 3 15 Oncorhynchus keta half-sib half-sib fish AN RR062b_b.G RR062b_b.N
./data_evolvability/RR062b_b/N2_n=98_P.csv Beacham 1988 RR062b_b N P 98 3 16 Oncorhynchus keta half-sib half-sib fish AN RR062b_b.P RR062b_b.N
./data_evolvability/RR073/B_n=68_G.csv Evans et al. 2015 RR073 B G 68 8 18 Poecilia reticulata half-sib half-sib fish BS RR073.G RR073.B
./data_evolvability/RR073/B_n=161_P.csv Evans et al. 2015 RR073 B P 161 8 17 Poecilia reticulata half-sib half-sib fish BS RR073.P RR073.B
./data_evolvability/RR073/S1_n=64_G.csv Evans et al. 2015 RR073 S G 64 8 20 Poecilia reticulata half-sib half-sib fish BS RR073.G RR073.S
./data_evolvability/RR073/S1_n=146_P.csv Evans et al. 2015 RR073 S P 146 8 19 Poecilia reticulata half-sib half-sib fish BS RR073.P RR073.S
./data_evolvability/RR076/B_n=69_G.csv Guan et al. 2016 RR076 B G 69 3 22 Scophthalmus maximus half-sib half-sib fish BS RR076.G RR076.B
./data_evolvability/RR076/B_n=2125_P.csv Guan et al. 2016 RR076 B P 2125 3 21 Scophthalmus maximus half-sib half-sib fish BS RR076.P RR076.B
./data_evolvability/RR076/S1_n=69_G.csv Guan et al. 2016 RR076 S G 69 3 24 Scophthalmus maximus half-sib half-sib fish BS RR076.G RR076.S
./data_evolvability/RR076/S1_n=2925_P.csv Guan et al. 2016 RR076 S P 2925 3 23 Scophthalmus maximus half-sib half-sib fish BS RR076.P RR076.S
./data_evolvability/RR087a/A_n=100_G.csv Sae-Lim et al. 2013 RR087a A G 100 3 25 Oncorhynchus mykiss full-sib full-sib fish AN RR087a.G RR087a.A
./data_evolvability/RR087a/A_n=2364_P.csv Sae-Lim et al. 2013 RR087a A P 2364 3 26 Oncorhynchus mykiss full-sib full-sib fish AN RR087a.P RR087a.A
./data_evolvability/RR087a/N1_n=100_G.csv Sae-Lim et al. 2013 RR087a N G 100 3 27 Oncorhynchus mykiss full-sib full-sib fish AN RR087a.G RR087a.N
./data_evolvability/RR087a/N1_n=1891_P.csv Sae-Lim et al. 2013 RR087a N P 1891 3 28 Oncorhynchus mykiss full-sib full-sib fish AN RR087a.P RR087a.N
./data_evolvability/RR087b/A_n=100_G.csv Sae-Lim et al. 2013 RR087b A G 100 3 29 Oncorhynchus mykiss full-sib full-sib fish AN RR087b.G RR087b.A
./data_evolvability/RR087b/A_n=2364_P.csv Sae-Lim et al. 2013 RR087b A P 2364 3 30 Oncorhynchus mykiss full-sib full-sib fish AN RR087b.P RR087b.A
./data_evolvability/RR087b/N2_n=100_G.csv Sae-Lim et al. 2013 RR087b N G 100 3 31 Oncorhynchus mykiss full-sib full-sib fish AN RR087b.G RR087b.N
./data_evolvability/RR087b/N2_n=2790_P.csv Sae-Lim et al. 2013 RR087b N P 2790 3 32 Oncorhynchus mykiss full-sib full-sib fish AN RR087b.P RR087b.N
./data_evolvability/RR087c/A_n=100_G.csv Sae-Lim et al. 2013 RR087c A G 100 3 33 Oncorhynchus mykiss full-sib full-sib fish AN RR087c.G RR087c.A
./data_evolvability/RR087c/A_n=2364_P.csv Sae-Lim et al. 2013 RR087c A P 2364 3 34 Oncorhynchus mykiss full-sib full-sib fish AN RR087c.P RR087c.A
./data_evolvability/RR087c/N3_n=100_G.csv Sae-Lim et al. 2013 RR087c N G 100 3 35 Oncorhynchus mykiss full-sib full-sib fish AN RR087c.G RR087c.N
./data_evolvability/RR087c/N3_n=1818_P.csv Sae-Lim et al. 2013 RR087c N P 1818 3 36 Oncorhynchus mykiss full-sib full-sib fish AN RR087c.P RR087c.N
./data_evolvability/sR097/B_n=26_G.csv Kasule 1991 sR097 B G 26 2 37 Dysdercus fasciatus half-sib half-sib insect BS sR097.G sR097.B
./data_evolvability/sR097/B_n=312_P.csv Kasule 1991 sR097 B P 312 2 38 Dysdercus fasciatus half-sib half-sib insect BS sR097.P sR097.B
./data_evolvability/sR097/S1_n=26_G.csv Kasule 1991 sR097 S G 26 2 39 Dysdercus fasciatus half-sib half-sib insect BS sR097.G sR097.S
./data_evolvability/sR097/S1_n=312_P.csv Kasule 1991 sR097 S P 312 2 40 Dysdercus fasciatus half-sib half-sib insect BS sR097.P sR097.S
./data_evolvability/sR110b/B_n=37_G.csv Blanckenhorn & Heyland 2004 sR110b B G 37 2 42 Scathophaga stercoraria half-sib half-sib insect BS sR110b.G sR110b.B
./data_evolvability/sR110b/B_n=128_P.csv Blanckenhorn & Heyland 2004 sR110b B P 128 2 41 Scathophaga stercoraria half-sib half-sib insect BS sR110b.P sR110b.B
./data_evolvability/sR110b/S1_n=27_G.csv Blanckenhorn & Heyland 2004 sR110b S G 27 2 43 Scathophaga stercoraria half-sib half-sib insect BS sR110b.G sR110b.S
./data_evolvability/sR110b/S1_n=98_P.csv Blanckenhorn & Heyland 2004 sR110b S P 98 2 44 Scathophaga stercoraria half-sib half-sib insect BS sR110b.P sR110b.S
./data_evolvability/sR110c/B_n=40_G.csv Blanckenhorn & Heyland 2004 sR110c B G 40 2 46 Scathophaga stercoraria half-sib half-sib insect BS sR110c.G sR110c.B
./data_evolvability/sR110c/B_n=171_P.csv Blanckenhorn & Heyland 2004 sR110c B P 171 2 45 Scathophaga stercoraria half-sib half-sib insect BS sR110c.P sR110c.B
./data_evolvability/sR110c/S1_n=33_G.csv Blanckenhorn & Heyland 2004 sR110c S G 33 2 48 Scathophaga stercoraria half-sib half-sib insect BS sR110c.G sR110c.S
./data_evolvability/sR110c/S1_n=129_P.csv Blanckenhorn & Heyland 2004 sR110c S P 129 2 47 Scathophaga stercoraria half-sib half-sib insect BS sR110c.P sR110c.S
./data_evolvability/sR115/A_n=350_G.csv Conner et al. 2003 sR115 A G 350 6 50 Raphanus raphanistrum half-sib half-sib plant AN sR115.G sR115.A
./data_evolvability/sR115/A_n=1133_P.csv Conner et al. 2003 sR115 A P 1133 6 49 Raphanus raphanistrum half-sib half-sib plant AN sR115.P sR115.A
./data_evolvability/sR115/N_n=593_G.csv Conner et al. 2003 sR115 N G 593 6 51 Raphanus raphanistrum half-sib half-sib plant AN sR115.G sR115.N
./data_evolvability/sR115/N_n=888_P.csv Conner et al. 2003 sR115 N P 888 6 52 Raphanus raphanistrum half-sib half-sib plant AN sR115.P sR115.N
./data_evolvability/sR116/A_n=366_G.csv Delcourt et al. 2010 sR116 A G 366 8 54 Drosophila serrata half-sib half-sib insect AN sR116.G sR116.A
./data_evolvability/sR116/A_n=3315_P.csv Delcourt et al. 2010 sR116 A P 3315 8 53 Drosophila serrata half-sib half-sib insect AN sR116.P sR116.A
./data_evolvability/sR116/N1_n=365_G.csv Delcourt et al. 2010 sR116 N G 365 8 56 Drosophila serrata half-sib half-sib insect AN sR116.G sR116.N
./data_evolvability/sR116/N1_n=3283_P.csv Delcourt et al. 2010 sR116 N P 3283 8 55 Drosophila serrata half-sib half-sib insect AN sR116.P sR116.N
./data_evolvability/sR145/B_n=52_G.csv Tucic et al. 1991 sR145 B G 52 8 58 Acanthoscelides obtectus full-sib full-sib insect BS sR145.G sR145.B
./data_evolvability/sR145/B_n=150_P.csv Tucic et al. 1991 sR145 B P 150 8 57 Acanthoscelides obtectus full-sib full-sib insect BS sR145.P sR145.B
./data_evolvability/sR145/S1_n=43_G.csv Tucic et al. 1991 sR145 S G 43 8 60 Acanthoscelides obtectus full-sib full-sib insect BS sR145.G sR145.S
./data_evolvability/sR145/S1_n=124_P.csv Tucic et al. 1991 sR145 S P 124 8 59 Acanthoscelides obtectus full-sib full-sib insect BS sR145.P sR145.S
./data_evolvability/WB001/B_n=30_G.csv Auld 2010 WB001 B G 30 3 62 Physa acuta full-sib full-sib mollusc BS WB001.G WB001.B
./data_evolvability/WB001/B_n=210_P.csv Auld 2010 WB001 B P 210 3 61 Physa acuta full-sib full-sib mollusc BS WB001.P WB001.B
./data_evolvability/WB001/S1_n=30_G.csv Auld 2010 WB001 S G 30 3 64 Physa acuta full-sib full-sib mollusc BS WB001.G WB001.S
./data_evolvability/WB001/S1_n=150_P.csv Auld 2010 WB001 S P 150 3 63 Physa acuta full-sib full-sib mollusc BS WB001.P WB001.S
./data_evolvability/WB001_b/B_n=30_G.csv Auld 2010 WB001_b B G 30 3 66 Physa acuta full-sib full-sib mollusc BS WB001_b.G WB001_b.B
./data_evolvability/WB001_b/B_n=210_P.csv Auld 2010 WB001_b B P 210 3 65 Physa acuta full-sib full-sib mollusc BS WB001_b.P WB001_b.B
./data_evolvability/WB001_b/S2_n=30_G.csv Auld 2010 WB001_b S G 30 3 68 Physa acuta full-sib full-sib mollusc BS WB001_b.G WB001_b.S
./data_evolvability/WB001_b/S2_n=210_P.csv Auld 2010 WB001_b S P 210 3 67 Physa acuta full-sib full-sib mollusc BS WB001_b.P WB001_b.S
./data_evolvability/WB001_c/B_n=30_G.csv Auld 2010 WB001_c B G 30 3 70 Physa acuta full-sib full-sib mollusc BS WB001_c.G WB001_c.B
./data_evolvability/WB001_c/B_n=210_P.csv Auld 2010 WB001_c B P 210 3 69 Physa acuta full-sib full-sib mollusc BS WB001_c.P WB001_c.B
./data_evolvability/WB001_c/S3_n=30_G.csv Auld 2010 WB001_c S G 30 3 72 Physa acuta full-sib full-sib mollusc BS WB001_c.G WB001_c.S
./data_evolvability/WB001_c/S3_n=150_P.csv Auld 2010 WB001_c S P 150 3 71 Physa acuta full-sib full-sib mollusc BS WB001_c.P WB001_c.S
./data_evolvability/WB002a/A_n=39_G.csv Begin & Roff 2001 WB002a A G 39 5 73 Gryllus pennsylvanicus full-sib full-sib insect AN WB002a.G WB002a.A
./data_evolvability/WB002a/A_n=705_P.csv Begin & Roff 2001 WB002a A P 705 5 74 Gryllus pennsylvanicus full-sib full-sib insect AN WB002a.P WB002a.A
./data_evolvability/WB002a/N1_n=39_G.csv Begin & Roff 2001 WB002a N G 39 5 75 Gryllus pennsylvanicus full-sib full-sib insect AN WB002a.G WB002a.N
./data_evolvability/WB002a/N1_n=505_P.csv Begin & Roff 2001 WB002a N P 505 5 76 Gryllus pennsylvanicus full-sib full-sib insect AN WB002a.P WB002a.N
./data_evolvability/WB003/A_n=56_G.csv Begin et al. 2004 WB003 A G 56 5 78 Gryllus firmus full-sib full-sib insect AN WB003.G WB003.A
./data_evolvability/WB003/A_n=533_P.csv Begin et al. 2004 WB003 A P 533 5 77 Gryllus firmus full-sib full-sib insect AN WB003.P WB003.A
./data_evolvability/WB003/N1_n=62_G.csv Begin et al. 2004 WB003 N G 62 5 79 Gryllus firmus full-sib full-sib insect AN WB003.G WB003.N
./data_evolvability/WB003/N1_n=862_P.csv Begin et al. 2004 WB003 N P 862 5 80 Gryllus firmus full-sib full-sib insect AN WB003.P WB003.N
./data_evolvability/WB003_2/A_n=56_G.csv Begin et al. 2004 WB003_2 A G 56 5 82 Gryllus firmus full-sib full-sib insect AN WB003_2.G WB003_2.A
./data_evolvability/WB003_2/A_n=533_P.csv Begin et al. 2004 WB003_2 A P 533 5 81 Gryllus firmus full-sib full-sib insect AN WB003_2.P WB003_2.A
./data_evolvability/WB003_2/N2_n=60_G.csv Begin et al. 2004 WB003_2 N G 60 5 83 Gryllus firmus full-sib full-sib insect AN WB003_2.G WB003_2.N
./data_evolvability/WB003_2/N2_n=715_P.csv Begin et al. 2004 WB003_2 N P 715 5 84 Gryllus firmus full-sib full-sib insect AN WB003_2.P WB003_2.N
./data_evolvability/WB003a/A_n=45_G.csv Begin et al. 2004 WB003a A G 45 5 86 Gryllus firmus full-sib full-sib insect AN WB003a.G WB003a.A
./data_evolvability/WB003a/A_n=262_P.csv Begin et al. 2004 WB003a A P 262 5 85 Gryllus firmus full-sib full-sib insect AN WB003a.P WB003a.A
./data_evolvability/WB003a/N1_n=41_G.csv Begin et al. 2004 WB003a N G 41 5 88 Gryllus firmus full-sib full-sib insect AN WB003a.G WB003a.N
./data_evolvability/WB003a/N1_n=315_P.csv Begin et al. 2004 WB003a N P 315 5 87 Gryllus firmus full-sib full-sib insect AN WB003a.P WB003a.N
./data_evolvability/WB003a_2/A_n=45_G.csv Begin et al. 2004 WB003a_2 A G 45 5 90 Gryllus firmus full-sib full-sib insect AN WB003a_2.G WB003a_2.A
./data_evolvability/WB003a_2/A_n=262_P.csv Begin et al. 2004 WB003a_2 A P 262 5 89 Gryllus firmus full-sib full-sib insect AN WB003a_2.P WB003a_2.A
./data_evolvability/WB003a_2/N2_n=32_G.csv Begin et al. 2004 WB003a_2 N G 32 5 92 Gryllus firmus full-sib full-sib insect AN WB003a_2.G WB003a_2.N
./data_evolvability/WB003a_2/N2_n=133_P.csv Begin et al. 2004 WB003a_2 N P 133 5 91 Gryllus firmus full-sib full-sib insect AN WB003a_2.P WB003a_2.N
./data_evolvability/WB007/B_n=132_G.csv Brock et al. 2010 WB007 B G 132 11 93 Brassica rapa RIL half-sib plant BS WB007.G WB007.B
./data_evolvability/WB007/B_n=601_P.csv Brock et al. 2010 WB007 B P 601 11 94 Brassica rapa RIL half-sib plant BS WB007.P WB007.B
./data_evolvability/WB007/S1_n=131_G.csv Brock et al. 2010 WB007 S G 131 11 95 Brassica rapa RIL half-sib plant BS WB007.G WB007.S
./data_evolvability/WB007/S1_n=524_P.csv Brock et al. 2010 WB007 S P 524 11 96 Brassica rapa RIL half-sib plant BS WB007.P WB007.S
./data_evolvability/WB010/A_n=141_G.csv Collins et al. 1999 WB010 A G 141 4 97 Achroia grisella half-sib half-sib insect AN WB010.G WB010.A
./data_evolvability/WB010/A_n=1540_P.csv Collins et al. 1999 WB010 A P 1540 4 98 Achroia grisella half-sib half-sib insect AN WB010.P WB010.A
./data_evolvability/WB010/N1_n=62_G.csv Collins et al. 1999 WB010 N G 62 4 100 Achroia grisella half-sib half-sib insect AN WB010.G WB010.N
./data_evolvability/WB010/N1_n=443_P.csv Collins et al. 1999 WB010 N P 443 4 99 Achroia grisella half-sib half-sib insect AN WB010.P WB010.N
./data_evolvability/WB011/B_n=404_G.csv Czesak & Fox 2003 WB011 B G 404 2 102 Stator limbatus half-sib half-sib insect BS WB011.G WB011.B
./data_evolvability/WB011/B_n=2667_P.csv Czesak & Fox 2003 WB011 B P 2667 2 101 Stator limbatus half-sib half-sib insect BS WB011.P WB011.B
./data_evolvability/WB011/S1_n=404_G.csv Czesak & Fox 2003 WB011 S G 404 2 104 Stator limbatus half-sib half-sib insect BS WB011.G WB011.S
./data_evolvability/WB011/S1_n=2667_P.csv Czesak & Fox 2003 WB011 S P 2667 2 103 Stator limbatus half-sib half-sib insect BS WB011.P WB011.S
./data_evolvability/WB014/A_n=1007_P.csv Delcourt et al. 2009 WB014 A P 1007 2 105 Drosophila serrata half-sib half-sib insect AN WB014.P WB014.A
./data_evolvability/WB014/A_n=319_G.csv Delcourt et al. 2009 WB014 A G 319 2 106 Drosophila serrata half-sib half-sib insect AN WB014.G WB014.A
./data_evolvability/WB014/N1_n=317_G.csv Delcourt et al. 2009 WB014 N G 317 2 107 Drosophila serrata half-sib half-sib insect AN WB014.G WB014.N
./data_evolvability/WB014/N1_n=998_P.csv Delcourt et al. 2009 WB014 N P 998 2 108 Drosophila serrata half-sib half-sib insect AN WB014.P WB014.N
./data_evolvability/WB016/B_n=27_G.csv Engqvist 2007 WB016 B G 27 2 110 Panorpa cognata full-sib full-sib insect BS WB016.G WB016.B
./data_evolvability/WB016/B_n=140_P.csv Engqvist 2007 WB016 B P 140 2 109 Panorpa cognata full-sib full-sib insect BS WB016.P WB016.B
./data_evolvability/WB016/S1_n=27_G.csv Engqvist 2007 WB016 S G 27 2 112 Panorpa cognata full-sib full-sib insect BS WB016.G WB016.S
./data_evolvability/WB016/S1_n=140_P.csv Engqvist 2007 WB016 S P 140 2 111 Panorpa cognata full-sib full-sib insect BS WB016.P WB016.S
./data_evolvability/WB022/A_n=40_G.csv Guntrip et al. 1997 WB022 A G 40 2 113 Callosobruchus maculatus half-sib half-sib insect AN WB022.G WB022.A
./data_evolvability/WB022/A_n=531_P.csv Guntrip et al. 1997 WB022 A P 531 2 114 Callosobruchus maculatus half-sib half-sib insect AN WB022.P WB022.A
./data_evolvability/WB022/N1_n=40_G.csv Guntrip et al. 1997 WB022 N G 40 2 115 Callosobruchus maculatus half-sib half-sib insect AN WB022.G WB022.N
./data_evolvability/WB022/N1_n=505_P.csv Guntrip et al. 1997 WB022 N P 505 2 116 Callosobruchus maculatus half-sib half-sib insect AN WB022.P WB022.N
./data_evolvability/WB025/A_n=65_G.csv Holloway et al. 1990 WB025 A G 65 4 117 Sitophilus oryzae half-sib half-sib insect AN WB025.G WB025.A
./data_evolvability/WB025/A_n=823_P.csv Holloway et al. 1990 WB025 A P 823 4 118 Sitophilus oryzae half-sib half-sib insect AN WB025.P WB025.A
./data_evolvability/WB025/N1_n=69_G.csv Holloway et al. 1990 WB025 N G 69 4 119 Sitophilus oryzae half-sib half-sib insect AN WB025.G WB025.N
./data_evolvability/WB025/N1_n=778_P.csv Holloway et al. 1990 WB025 N P 778 4 120 Sitophilus oryzae half-sib half-sib insect AN WB025.P WB025.N
./data_evolvability/WB028/B_n=60_G.csv Klause & Morin 2001 WB028 B G 60 2 121 Priophorus pallipes full-sib full-sib insect BS WB028.G WB028.B
./data_evolvability/WB028/B_n=600_P.csv Klause & Morin 2001 WB028 B P 600 2 122 Priophorus pallipes full-sib full-sib insect BS WB028.P WB028.B
./data_evolvability/WB028/S1_n=60_G.csv Klause & Morin 2001 WB028 S G 60 2 123 Priophorus pallipes full-sib full-sib insect BS WB028.G WB028.S
./data_evolvability/WB028/S1_n=600_P.csv Klause & Morin 2001 WB028 S P 600 2 124 Priophorus pallipes full-sib full-sib insect BS WB028.P WB028.S
./data_evolvability/WB029/B_n=63_G.csv King et al. 2011 WB029 B G 63 3 126 Gryllus firmus half-sib half-sib insect BS WB029.G WB029.B
./data_evolvability/WB029/B_n=1662_P.csv King et al. 2011 WB029 B P 1662 3 125 Gryllus firmus half-sib half-sib insect BS WB029.P WB029.B
./data_evolvability/WB029/S1_n=63_G.csv King et al. 2011 WB029 S G 63 3 128 Gryllus firmus half-sib half-sib insect BS WB029.G WB029.S
./data_evolvability/WB029/S1_n=1728_P.csv King et al. 2011 WB029 S P 1728 3 127 Gryllus firmus half-sib half-sib insect BS WB029.P WB029.S
./data_evolvability/WB033/A_n=60_G.csv Lau et al. 2014 WB033 A G 60 4 130 Arabidopsis thaliana RIL half-sib plant AN WB033.G WB033.A
./data_evolvability/WB033/A_n=526_P.csv Lau et al. 2014 WB033 A P 526 4 129 Arabidopsis thaliana RIL half-sib plant AN WB033.P WB033.A
./data_evolvability/WB033/N1_n=60_G.csv Lau et al. 2014 WB033 N G 60 4 132 Arabidopsis thaliana RIL half-sib plant AN WB033.G WB033.N
./data_evolvability/WB033/N1_n=500_P.csv Lau et al. 2014 WB033 N P 500 4 131 Arabidopsis thaliana RIL half-sib plant AN WB033.P WB033.N
./data_evolvability/WB033a/A_n=60_G.csv Lau et al. 2014 WB033a A G 60 4 134 Arabidopsis thaliana RIL half-sib plant AN WB033a.G WB033a.A
./data_evolvability/WB033a/A_n=459_P.csv Lau et al. 2014 WB033a A P 459 4 133 Arabidopsis thaliana RIL half-sib plant AN WB033a.P WB033a.A
./data_evolvability/WB033a/N1_n=60_G.csv Lau et al. 2014 WB033a N G 60 4 136 Arabidopsis thaliana RIL half-sib plant AN WB033a.G WB033a.N
./data_evolvability/WB033a/N1_n=470_P.csv Lau et al. 2014 WB033a N P 470 4 135 Arabidopsis thaliana RIL half-sib plant AN WB033a.P WB033a.N
./data_evolvability/WB035/B_n=432_G.csv Messina et al. 2013 WB035 B G 432 3 138 Callosobruchus maculatus half-sib half-sib insect BS WB035.G WB035.B
./data_evolvability/WB035/B_n=2201_P.csv Messina et al. 2013 WB035 B P 2201 3 137 Callosobruchus maculatus half-sib half-sib insect BS WB035.P WB035.B
./data_evolvability/WB035/S1_n=432_G.csv Messina et al. 2013 WB035 S G 432 3 140 Callosobruchus maculatus half-sib half-sib insect BS WB035.G WB035.S
./data_evolvability/WB035/S1_n=2207_P.csv Messina et al. 2013 WB035 S P 2207 3 139 Callosobruchus maculatus half-sib half-sib insect BS WB035.P WB035.S
./data_evolvability/WB036/A_n=22_G.csv Paccard et al. 2013 WB036 A G 22 7 141 Arabidopsis lyrata full-sib full-sib plant AN WB036.G WB036.A
./data_evolvability/WB036/A_n=58_P.csv Paccard et al. 2013 WB036 A P 58 7 142 Arabidopsis lyrata full-sib full-sib plant AN WB036.P WB036.A
./data_evolvability/WB036/N1_n=22_G.csv Paccard et al. 2013 WB036 N G 22 7 143 Arabidopsis lyrata full-sib full-sib plant AN WB036.G WB036.N
./data_evolvability/WB036/N1_n=63_P.csv Paccard et al. 2013 WB036 N P 63 7 144 Arabidopsis lyrata full-sib full-sib plant AN WB036.P WB036.N
./data_evolvability/WB036A/A_n=22_G.csv Paccard et al. 2013 WB036A A G 22 7 145 Arabidopsis lyrata full-sib full-sib plant AN WB036A.G WB036A.A
./data_evolvability/WB036A/A_n=60_P.csv Paccard et al. 2013 WB036A A P 60 7 146 Arabidopsis lyrata full-sib full-sib plant AN WB036A.P WB036A.A
./data_evolvability/WB036A/N1_n=22_G.csv Paccard et al. 2013 WB036A N G 22 7 147 Arabidopsis lyrata full-sib full-sib plant AN WB036A.G WB036A.N
./data_evolvability/WB036A/N1_n=63_P.csv Paccard et al. 2013 WB036A N P 63 7 148 Arabidopsis lyrata full-sib full-sib plant AN WB036A.P WB036A.N
./data_evolvability/WB038/B_n=135_G.csv Rauter&Moore 2002 WB038 B G 135 4 149 Nicrophorus pustulatus half-sib half-sib insect BS WB038.G WB038.B
./data_evolvability/WB038/B_n=639_P.csv Rauter&Moore 2002 WB038 B P 639 4 150 Nicrophorus pustulatus half-sib half-sib insect BS WB038.P WB038.B
./data_evolvability/WB038/S1_n=135_G.csv Rauter&Moore 2002 WB038 S G 135 4 151 Nicrophorus pustulatus half-sib half-sib insect BS WB038.G WB038.S
./data_evolvability/WB038/S1_n=636_P.csv Rauter&Moore 2002 WB038 S P 636 4 152 Nicrophorus pustulatus half-sib half-sib insect BS WB038.P WB038.S
./data_evolvability/WB039/B_n=21_G.csv Relyea 2005 WB039 B G 21 9 153 Rana sylvatica half-sib half-sib amphibian BS WB039.G WB039.B
./data_evolvability/WB039/B_n=2700_P.csv Relyea 2005 WB039 B P 2700 9 154 Rana sylvatica half-sib half-sib amphibian BS WB039.P WB039.B
./data_evolvability/WB039/S1_n=21_G.csv Relyea 2005 WB039 S G 21 9 155 Rana sylvatica half-sib half-sib amphibian BS WB039.G WB039.S
./data_evolvability/WB039/S1_n=2700_P.csv Relyea 2005 WB039 S P 2700 9 156 Rana sylvatica half-sib half-sib amphibian BS WB039.P WB039.S
./data_evolvability/WB044/B_n=26_G.csv Sherrard et al. 2009 WB044 B G 26 10 158 Avena barbata RIL half-sib plant BS WB044.G WB044.B
./data_evolvability/WB044/B_n=112_P.csv Sherrard et al. 2009 WB044 B P 112 10 157 Avena barbata RIL half-sib plant BS WB044.P WB044.B
./data_evolvability/WB044/S1_n=26_G.csv Sherrard et al. 2009 WB044 S G 26 10 160 Avena barbata RIL half-sib plant BS WB044.G WB044.S
./data_evolvability/WB044/S1_n=112_P.csv Sherrard et al. 2009 WB044 S P 112 10 159 Avena barbata RIL half-sib plant BS WB044.P WB044.S
./data_evolvability/WB045/A_n=39_G.csv Simons & Roff 1996 WB045 A G 39 7 162 Gryllus pennsylvanicus full-sib full-sib insect AN WB045.G WB045.A
./data_evolvability/WB045/A_n=380_P.csv Simons & Roff 1996 WB045 A P 380 7 161 Gryllus pennsylvanicus full-sib full-sib insect AN WB045.P WB045.A
./data_evolvability/WB045/N1_n=39_G.csv Simons & Roff 1996 WB045 N G 39 7 164 Gryllus pennsylvanicus full-sib full-sib insect AN WB045.G WB045.N
./data_evolvability/WB045/N1_n=362_P.csv Simons & Roff 1996 WB045 N P 362 7 163 Gryllus pennsylvanicus full-sib full-sib insect AN WB045.P WB045.N
./data_evolvability/WB045a/A_n=39_G.csv Simons & Roff 1996 WB045a A G 39 6 165 Gryllus pennsylvanicus full-sib full-sib insect AN WB045a.G WB045a.A
./data_evolvability/WB045a/A_n=391_P.csv Simons & Roff 1996 WB045a A P 391 6 166 Gryllus pennsylvanicus full-sib full-sib insect AN WB045a.P WB045a.A
./data_evolvability/WB045a/N1_n=39_G.csv Simons & Roff 1996 WB045a N G 39 6 168 Gryllus pennsylvanicus full-sib full-sib insect AN WB045a.G WB045a.N
./data_evolvability/WB045a/N1_n=371_P.csv Simons & Roff 1996 WB045a N P 371 6 167 Gryllus pennsylvanicus full-sib full-sib insect AN WB045a.P WB045a.N
./data_evolvability/WB046/B_n=50_G.csv Simonsen & Stinchcombe 2001 WB046 B G 50 3 169 Ipomoea hederacea RIL half-sib plant BS WB046.G WB046.B
./data_evolvability/WB046/B_n=600_P.csv Simonsen & Stinchcombe 2001 WB046 B P 600 3 170 Ipomoea hederacea RIL half-sib plant BS WB046.P WB046.B
./data_evolvability/WB046/S1_n=50_G.csv Simonsen & Stinchcombe 2001 WB046 S G 50 3 171 Ipomoea hederacea RIL half-sib plant BS WB046.G WB046.S
./data_evolvability/WB046/S1_n=600_P.csv Simonsen & Stinchcombe 2001 WB046 S P 600 3 172 Ipomoea hederacea RIL half-sib plant BS WB046.P WB046.S
./data_evolvability/WB053/A_n=24_G.csv Via & Conner 1995 WB053 A G 24 2 173 Tribolium castaneum half-sib half-sib insect AN WB053.G WB053.A
./data_evolvability/WB053/A_n=480_P.csv Via & Conner 1995 WB053 A P 480 2 174 Tribolium castaneum half-sib half-sib insect AN WB053.P WB053.A
./data_evolvability/WB053/N1_n=24_G.csv Via & Conner 1995 WB053 N G 24 2 175 Tribolium castaneum half-sib half-sib insect AN WB053.G WB053.N
./data_evolvability/WB053/N1_n=480_P.csv Via & Conner 1995 WB053 N P 480 2 176 Tribolium castaneum half-sib half-sib insect AN WB053.P WB053.N
./data_evolvability/WB053_3/A_n=24_G.csv Via & Conner 1995 WB053_3 A G 24 2 177 Tribolium castaneum half-sib half-sib insect AN WB053_3.G WB053_3.A
./data_evolvability/WB053_3/A_n=480_P.csv Via & Conner 1995 WB053_3 A P 480 2 178 Tribolium castaneum half-sib half-sib insect AN WB053_3.P WB053_3.A
./data_evolvability/WB053_3/N3_n=24_G.csv Via & Conner 1995 WB053_3 N G 24 2 179 Tribolium castaneum half-sib half-sib insect AN WB053_3.G WB053_3.N
./data_evolvability/WB053_3/N3_n=480_P.csv Via & Conner 1995 WB053_3 N P 480 2 180 Tribolium castaneum half-sib half-sib insect AN WB053_3.P WB053_3.N
./data_evolvability/WB053b/A_n=80_G.csv Via & Conner 1995 WB053b A G 80 2 182 Tribolium castaneum half-sib half-sib insect AN WB053b.G WB053b.A
./data_evolvability/WB053b/A_n=480_P.csv Via & Conner 1995 WB053b A P 480 2 181 Tribolium castaneum half-sib half-sib insect AN WB053b.P WB053b.A
./data_evolvability/WB053b/N1_n=80_G.csv Via & Conner 1995 WB053b N G 80 2 184 Tribolium castaneum half-sib half-sib insect AN WB053b.G WB053b.N
./data_evolvability/WB053b/N1_n=480_P.csv Via & Conner 1995 WB053b N P 480 2 183 Tribolium castaneum half-sib half-sib insect AN WB053b.P WB053b.N
./data_evolvability/WB053b_3/A_n=80_G.csv Via & Conner 1995 WB053b_3 A G 80 2 186 Tribolium castaneum half-sib half-sib insect AN WB053b_3.G WB053b_3.A
./data_evolvability/WB053b_3/A_n=480_P.csv Via & Conner 1995 WB053b_3 A P 480 2 185 Tribolium castaneum half-sib half-sib insect AN WB053b_3.P WB053b_3.A
./data_evolvability/WB053b_3/N3_n=80_G.csv Via & Conner 1995 WB053b_3 N G 80 2 188 Tribolium castaneum half-sib half-sib insect AN WB053b_3.G WB053b_3.N
./data_evolvability/WB053b_3/N3_n=480_P.csv Via & Conner 1995 WB053b_3 N P 480 2 187 Tribolium castaneum half-sib half-sib insect AN WB053b_3.P WB053b_3.N
./data_evolvability/WB056/A_n=64_G.csv Service & Rose 1985 WB056 A G 64 2 189 Drosophila melanogaster half-sib half-sib insect AN WB056.G WB056.A
./data_evolvability/WB056/A_n=758_P.csv Service & Rose 1985 WB056 A P 758 2 190 Drosophila melanogaster half-sib half-sib insect AN WB056.P WB056.A
./data_evolvability/WB056/N1_n=64_G.csv Service & Rose 1985 WB056 N G 64 2 191 Drosophila melanogaster half-sib half-sib insect AN WB056.G WB056.N
./data_evolvability/WB056/N1_n=758_P.csv Service & Rose 1985 WB056 N P 758 2 192 Drosophila melanogaster half-sib half-sib insect AN WB056.P WB056.N
./data_evolvability/WB057/A_n=22_G.csv Grill et al. 1997 WB057 A G 22 4 193 Harmonia axyridis full-sib full-sib insect AN WB057.G WB057.A
./data_evolvability/WB057/A_n=277_P.csv Grill et al. 1997 WB057 A P 277 4 194 Harmonia axyridis full-sib full-sib insect AN WB057.P WB057.A
./data_evolvability/WB057/N1_n=26_G.csv Grill et al. 1997 WB057 N G 26 4 195 Harmonia axyridis full-sib full-sib insect AN WB057.G WB057.N
./data_evolvability/WB057/N1_n=277_P.csv Grill et al. 1997 WB057 N P 277 4 196 Harmonia axyridis full-sib full-sib insect AN WB057.P WB057.N

We can then import the matrices and descriptive data.

         mats_evolve <- import_mat(matrices_evolve)
        dscrp_evolve <- import_mat(descripDat_evolve)

Scale the data to prevent disproportional large effects of the traits with large values:

  standard <- list()
  for (i in 1:length(dscrp_evolve)){
    # standardize descriptive data
    # calculate the average mean for each trait
    dscrp_evolve[[i]]$stanfac <- ave(dscrp_evolve[[i]]$mean,dscrp_evolve[[i]]$trait)
    #standardize mean, sd, se, Va
    dscrp_evolve[[i]]$mean <- dscrp_evolve[[i]]$mean/dscrp_evolve[[i]]$stanfac
    dscrp_evolve[[i]]$sd <- dscrp_evolve[[i]]$sd/dscrp_evolve[[i]]$stanfac
    dscrp_evolve[[i]]$se <- dscrp_evolve[[i]]$se/dscrp_evolve[[i]]$stanfac
    dscrp_evolve[[i]]$Va <- dscrp_evolve[[i]]$Va/(dscrp_evolve[[i]]$stanfac^2)
    # save the standardization factor in to another list
    standard[[i]] <- data.frame(trait=unique(dscrp_evolve[[i]]$trait),stanfac=dscrp_evolve[[i]]$stanfac[which(!duplicated(dscrp_evolve[[i]]$trait))]) 
    # remove spaces and dashes from the trait names to be consistent with the names in the covariance-matrices.
    standard[[i]]$trait <- gsub(" ","\\.", as.character(standard[[i]]$trait))
    standard[[i]]$trait <- gsub("-","\\.", standard[[i]]$trait)
    rownames(standard[[i]]) <- standard[[i]]$trait
    }
  
  id_d <- sapply(names(dscrp_evolve),function(x){strsplit(x,"/")[[1]][3]})
  id_m <- sapply(names(mats_evolve),function(x){strsplit(x,"/")[[1]][3]})
  for (i in 1:length(mats_evolve)){
    # select the right standardization factors
    j <- which(id_d == id_m[[i]])
    std <- standard[[j]]
    # remove traits not represented in the covar-matrix
    std <- std[which(std$trait %in% names(mats_evolve[[i]])),]
    mats_evolve[[i]] <- mats_evolve[[i]][which(names(mats_evolve[[i]]) %in% std$trait),]
    mats_evolve[[i]] <- mats_evolve[[i]][,which(names(mats_evolve[[i]]) %in% std$trait)]
    # order traits as in the covar-matrix
    std <- std[names(mats_evolve[[i]]),]
    # check whether the standardization factors line up
    if(sum(names(mats_evolve[[i]]) != std$trait)>0){stop("There is a inconsistency between trait names of a matrix and a description file")}
    mats_evolve[[i]] <- mats_evolve[[i]]/(std$stanfac %*% t(std$stanfac))
    }
  
  length(mats_evolve)
## [1] 196
  length(dscrp_evolve)
## [1] 49

2.2 – Descriptors of the Studies

Using these matrices we can also get some descriptor information about the studies

           env_mats_evolve <- extract_matrix_env(names(mats_evolve), type = "full")
         names_mats_evolve <- extract_matrix_studynames(names(mats_evolve))
          type_mats_evolve <- extract_matrix_type(names(mats_evolve))

        means_env <- plyr::ldply(lapply(dscrp_evolve, function(x) x[,c("env","trait","mean")]))
        
        # make trait names consistent with the row and column names
        means_env$trait <- gsub(" ","\\.",means_env$trait)
        means_env$trait <- gsub("_","\\.",means_env$trait)
        
    means_env$.id <- gsub("/DecsrpData.csv", "", means_env$.id)
    means_env$.id <- gsub("./data/", "", means_env$.id)
  means_env$env_r <- gsub("[1-9]", "", as.character(means_env$env))
  means_env$match <- paste0(means_env$.id, "_", means_env$env)

Now that we have extracted the mean vectors for all traits within each study and environment we can extract the G matrices for novel and non-novel environments.

           G_matrices <- mats_evolve[ grep("G", names(mats_evolve)) ]
        G_matrices_NS <- G_matrices[ grep("/[NS]", names(G_matrices)) ]
        G_matrices_AB <- G_matrices[ grep("/[AB]", names(G_matrices)) ]

        n_G_NS <- match_n(names(G_matrices_NS))
        n_G_AB <- match_n(names(G_matrices_AB))

We can check that everything matches, so we can be sure that nothing has been disorganized using some simple tests

    # First, we can subset just the G matrices from the metadata file
        tst <- meta_evolve %>% filter(type == "G")

    # Extract n from the G matrices
        n_tst <- match_n(names(G_matrices))

    # We can test that the sample sizes match what is actually in the file names
        all(tst$n == n_tst)  # TRUE, so we are good here.
## [1] FALSE
    # We can also make sure that the col matches with the G-matrices
        all(tst$col == match(tst$file, names(mats_evolve))) # TRUE, so good.
## [1] FALSE
    # Check that the sample sizes match with the matrices for each matrix. That way we can be absolutely sure that N matches the matrix
        n_NS_tst <- match_n(names(G_matrices_NS))
        n_AB_tst <- match_n(names(G_matrices_AB))

        all(n_G_NS == n_NS_tst) # TRUE
## [1] TRUE
        all(n_G_AB == n_AB_tst) # TRUE
## [1] TRUE

2.3 – Simulations for Evolvability Metrics

We now have 5000 simulated G matrices with the sample sizes (number of families) used to estimate these matrices. Given that we don’t have the P matrices for these data, we can’t simulate changes in the means using the same approaches (but see below). For now, we will use the mean differences in traits across environments to superimpose the plasticity vector on each G matrix. First, we need the centroid differences in means between the two environments:

    # To make things easy to compare, just split the data
    means_env$study <- gsub("./data_evolvability/", "", means_env$.id)
    split_mean_stdy <- base::split(means_env, means_env$study)
    split_mean_stdy <- lapply(split_mean_stdy, function(x) arrange(x, env))

    # Get the plasticity response vector between the non-novel / benign and novel / stressful.
    cen.diff <- lapply(split_mean_stdy, function(x)centroid.diff(x[x$env_r == "A" | x$env_r == "B",]$mean, x[x$env_r == "N" | x$env_r == "S",]$mean)$vector)
    
    # Get corresponding trait names
    cen.diff.names <- lapply(split_mean_stdy, function(x){as.character(x[x$env_r == "A" | x$env_r == "B",]$trait)})

Locate cases in which the number of traits do not match up

    ids <- unique(c(which(unlist(lapply(cen.diff,length))!=unlist(lapply(G_matrices_AB,function(x){dim(x)[1]}))),which(unlist(lapply(cen.diff,length))!=unlist(lapply(G_matrices_NS,function(x){dim(x)[1]})))))

    for(i in ids){
      # Identify the trait which all have in common
      intsect <- intersect(cen.diff.names[[i]],intersect(names(G_matrices_NS[[i]]),names(G_matrices_AB[[i]])))
      # select only the traits in common
      cen.diff[[i]] <- cen.diff[[i]][which(cen.diff.names[[i]]%in%intsect)]
      G_matrices_AB[[i]] <- G_matrices_AB[[i]][which(names(G_matrices_AB[[i]])%in%intsect),]
      G_matrices_AB[[i]] <- G_matrices_AB[[i]][,which(names(G_matrices_AB[[i]])%in%intsect)]
      G_matrices_NS[[i]] <- G_matrices_NS[[i]][which(names(G_matrices_NS[[i]])%in%intsect),]
      G_matrices_NS[[i]] <- G_matrices_NS[[i]][,which(names(G_matrices_NS[[i]])%in%intsect)]  
    }

Now that we have checked things and we have the novel (i.e., NS) and non-novel (i.e., AB) matrices extracted we can do some simulations with the matrices

    # Simulate G matrices for environments
        sims = 5000

        G_sim_NS <- mapply(function(x,y) sim_G(x, y, sims = sims, cor = FALSE, matrix_fix = "bend"), 
                            x = G_matrices_NS, 
                            y = split(n_G_NS, 1:length(n_G_NS)))
        G_sim_AB <- mapply(function(x,y) sim_G(x, y, sims = sims, cor = FALSE, matrix_fix = "bend"), 
                            x = G_matrices_AB, 
                            y = split(n_G_AB, 1:length(n_G_AB)))

Plasticity vectors are based on the multivariate mean centroid differences. We can now project these vectors on the G matrices in the non-novel and novel environment to understand how plastic responses brought about by the environment lead to greater, lower or no change in evolvability. We note that the trait means themselves also have sampling variances associated with them, however, we are assuming that these are fixed given we don’t have a P matrix to simulate the mean response vector for all studies. However, given that we are incorporating sampling error associated with G this will drive weighted differences among studies. Adding any additional sampling variance resulting from P vectors is likely not to change results because this will just add proportionally more sampling variance across studies in many cases.

    # Here we can get the amount of genetic variance along the plastic vector
        evolv_max_G_NS <- evolveability_along_plastic_vector(G_sim_NS, cen.diff, compare = "max")

        evolv_max_G_AB <- evolveability_along_plastic_vector(G_sim_AB, cen.diff, compare = "max")

    # Here we test whether the amount of genetic variance along the plastic vector is higher than expected
        evolv_max_G_NS_test <- evolveability_along_plastic_vector(G_sim_NS, cen.diff, compare = "maxcorr")

        evolv_max_G_AB_test <- evolveability_along_plastic_vector(G_sim_AB, cen.diff, compare = "maxcorr")
        
    # Here we can get the angle between the plastic vector and gmax
        evol_angle_NS <- evolveability_angle_along_plastic_vector(G_sim_NS, cen.diff)
        
        evol_angle_AB <- evolveability_angle_along_plastic_vector(G_sim_AB, cen.diff)

2.4 – Plotting Simulations

Now that we have all the simulations finished we can explore the distribution of effects across the studies

    plot_effect_dist(evolv_max_G_NS)

Figure 2.4.1 – Distribution of effects for the proportion of genetic variation along plastic response vector for \(G_{novel}\).

    plot_effect_dist(evolv_max_G_AB)

Figure 2.4.2 – Distribution of effects for the proportion of genetic variation along plastic response vector for \(G_{non-novel}\).

    plot_effect_dist(evol_angle_NS)

Figure 2.4.3 – Distribution of effects for the angle between \(gmax_{novel}\) and the plastic response vector for \(G_{novel}\).

    plot_effect_dist(evol_angle_AB)

Figure 2.4.4 – Distribution of effects for the angle between \(gmax_{nonnovel}\) and the plastic response vector for \(G_{non-novel}\).

2.5 – Generating the Data for Meta-analysis

We now want to condense down the simulated distributions and sampling variance so that we can more easily meta-analyse these data. Of course, one could technically use the entire distribution, but computationally this will make the models take a long time to run. So, we will keep it simple. We will first generate the data for \(\pi_{e}\).

                effect_G_NS <- create_data_evolve(evolv_max_G_NS, meta_evolve, env = c("N", "S"))

                effect_G_AB <- create_data_evolve(evolv_max_G_AB, meta_evolve, env = c("A", "B"))

            evolve_max_data <- rbind(effect_G_NS, effect_G_AB)
        evolve_max_data$err <- 1:nrow(evolve_max_data)
evolve_max_data$environment <- ifelse(evolve_max_data$env == "B" | evolve_max_data$env == "A", "Non-novel", "Novel")
         evolve_max_data$stdy <- substring(evolve_max_data$study,1, 5)
       
        write.csv(evolve_max_data, file = "./data_evolvability/evolve_max_data_standardised.csv",row.names=FALSE)

Next for will generate the data for \(\pi_{e-0}\).

                effect_G_NS_test <- create_data_evolve(evolv_max_G_NS_test, meta_evolve, env = c("N", "S"))

                effect_G_AB_test <- create_data_evolve(evolv_max_G_AB_test, meta_evolve, env = c("A", "B"))

            evolve_max_data_test <- rbind(effect_G_NS_test, effect_G_AB_test)
        evolve_max_data_test$err <- 1:nrow(evolve_max_data_test)
evolve_max_data_test$environment <- ifelse(evolve_max_data_test$env == "B" | evolve_max_data_test$env == "A", "Non-novel", "Novel")
         evolve_max_data_test$stdy <- substring(evolve_max_data_test$study,1, 5)
       
        write.csv(evolve_max_data_test, file = "./data_evolvability/evolve_max_data_test.csv",row.names=FALSE)

Now that we have data for \(pi_{e}\), we can generate \(\theta_{e}\) in the same way:

                     angle_G_NS <- create_data_evolve(evol_angle_NS, meta_evolve, env = c("N", "S"))

                     angle_G_AB <- create_data_evolve(evol_angle_AB, meta_evolve, env = c("A", "B"))

              evolve_angle_data <- rbind(angle_G_NS, angle_G_AB)
          evolve_angle_data$err <- 1:nrow(evolve_angle_data)
  evolve_angle_data$environment <- ifelse(evolve_angle_data$env == "B" | evolve_angle_data$env == "A", "Non-novel", "Novel")
       evolve_angle_data$stdy <- substring(evolve_angle_data$study,1, 5)


        write.csv(evolve_angle_data, file = "./data_evolvability/evolve_angle_data_standardised.csv", row.names=FALSE)

3.0 – Meta-analysis 1: General patterns of G and P within and across environments

3.1 – Load Data

    # Load libraries
        library(metafor)
        library(knitr)
        library(rmarkdown)

    # Load data
                SMD_data <- read.csv("./data_GP_compare/SMD_data_standardised.csv")
              angle_data <- read.csv("./data_GP_compare/angle_across_data_standardised.csv")
       angle_data_within <- read.csv("./data_GP_compare/angle_within_data_standardised.csv")
    var_along_logit_data <- read.csv("./data_GP_compare/var_along_logit_data_standardised.csv")

3.2 – Meta-analysis of SDV

First, we can have a look at the dataset more carefully before running the analysis.

    kable(SMD_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")
X study type ES S_var authors year avg_n_type trait_num genus species design group env.comp species_full stdy err wi
1 RR073 G -1.1114684 0.0137012 Evans et al. 2015 66.0 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 1 72.98608
2 RR073 P -1.1761978 0.0036433 Evans et al. 2015 153.5 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 2 274.47870
3 RR076 G -1.3864048 0.0141709 Guan et al. 2016 69.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 3 70.56698
4 RR076 P -0.4165704 0.0010826 Guan et al. 2016 2525.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 4 923.73149
5 RR087a G -0.0463060 0.0266880 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 5 37.47005
6 RR087a P -0.1584050 0.0011757 Sae-Lim et al. 2013 2127.5 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 6 850.52360
7 RR087b G -0.0811399 0.0251500 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 7 39.76148
8 RR087b P 0.0924924 0.0009454 Sae-Lim et al. 2013 2577.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 8 1057.73026
9 RR087c G -0.2067433 0.0217486 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 9 45.97998
10 RR087c P -0.0052831 0.0010939 Sae-Lim et al. 2013 2091.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 10 914.14934
11 sR097 G -1.0757071 0.0647484 Kasule 1991 26.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 11 15.44440
12 sR097 P 0.5490437 0.0073155 Kasule 1991 312.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 12 136.69640
13 sR110b G 0.2457675 0.0971304 Blanckenhorn & Heyland 2004 32.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 13 10.29544
14 sR110b P 0.1302658 0.0254568 Blanckenhorn & Heyland 2004 113.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 14 39.28231
15 sR110c G -0.1644129 0.0885876 Blanckenhorn & Heyland 2004 36.5 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 15 11.28826
16 sR110c P -0.2031563 0.0200968 Blanckenhorn & Heyland 2004 150.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 16 49.75914
17 sR115 G -0.9262440 0.0020576 Conner et al. 2003 471.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 17 485.99435
18 sR115 P 0.4498877 0.0017084 Conner et al. 2003 1010.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 18 585.34934
19 sR116 G -0.3003715 0.0069157 Delcourt et al. 2010 365.5 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 19 144.59791
20 sR116 P -0.0686755 0.0009931 Delcourt et al. 2010 3299.0 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 20 1006.93838
21 sR145 G 0.1496778 0.0308049 Tucic et al. 1991 47.5 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 21 32.46242
22 sR145 P 0.0567341 0.0116437 Tucic et al. 1991 137.0 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 22 85.88331
23 WB002a G -0.0332739 0.0675742 Begin & Roff 2001 39.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 23 14.79854
24 WB002a P -0.2989204 0.0033918 Begin & Roff 2001 605.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 24 294.82803
25 WB003 G 0.8505977 0.0343841 Begin et al. 2004 59.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 25 29.08320
26 WB003 P 0.4246671 0.0038752 Begin et al. 2004 697.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 26 258.05423
27 WB003_2 G 0.3356583 0.0435435 Begin et al. 2004 58.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 27 22.96552
28 WB003_2 P 0.3286113 0.0043748 Begin et al. 2004 624.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 28 228.58439
29 WB003a G -0.0378720 0.0682350 Begin et al. 2004 43.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 29 14.65523
30 WB003a P 0.3435527 0.0093449 Begin et al. 2004 288.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 30 107.01026
31 WB003a_2 G 0.7981247 0.0636134 Begin et al. 2004 38.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 31 15.71997
32 WB003a_2 P 0.4089151 0.0164377 Begin et al. 2004 197.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 32 60.83565
33 WB007 G -1.0716421 0.0068572 Brock et al. 2010 131.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 33 145.83305
34 WB007 P -1.1222659 0.0011237 Brock et al. 2010 562.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 34 889.90631
35 WB010 G 0.1409895 0.0374687 Collins et al. 1999 101.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 35 26.68894
36 WB010 P 0.0187338 0.0042686 Collins et al. 1999 991.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 36 234.26914
37 WB011 G 0.2529621 0.0095757 Czesak & Fox 2003 404.0 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 37 104.43115
38 WB011 P 0.0000158 0.0010038 Czesak & Fox 2003 2667.0 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 38 996.20709
39 WB028 G 1.3743055 0.0162181 Klause & Morin 2001 60.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 39 61.65935
40 WB028 P 1.4176906 0.0012685 Klause & Morin 2001 600.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 40 788.30381
41 WB029 G -0.2907661 0.0601083 King et al. 2011 63.0 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 41 16.63664
42 WB029 P -1.0680634 0.0008849 King et al. 2011 1695.0 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 42 1130.05446
43 WB033 G 0.0528833 0.0663470 Lau et al. 2014 60.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 43 15.07227
44 WB033 P 0.5024050 0.0069455 Lau et al. 2014 513.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 44 143.97718
45 WB033a G 0.9189091 0.0414328 Lau et al. 2014 60.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 45 24.13545
46 WB033a P -0.0210253 0.0089497 Lau et al. 2014 464.5 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 46 111.73533
47 WB035 G -0.0573465 0.0064297 Messina et al. 2013 432.0 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 47 155.52774
48 WB035 P -1.2368324 0.0004308 Messina et al. 2013 2204.0 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 48 2321.03330
49 WB038 G -0.0589241 0.0177160 Rauter&Moore 2002 135.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 49 56.44611
50 WB038 P -0.3677845 0.0032270 Rauter&Moore 2002 639.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 50 309.88875
51 WB045 G 1.6727335 0.0085023 Simons & Roff 1996 39.0 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 51 117.61495
52 WB045 P 0.2529547 0.0076488 Simons & Roff 1996 371.0 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 52 130.73914
53 WB045a G 0.5721094 0.0636117 Simons & Roff 1996 39.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 53 15.72038
54 WB045a P 0.1596246 0.0067022 Simons & Roff 1996 381.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 54 149.20430
55 WB053 G -1.3688710 0.0476525 Via & Conner 1995 24.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 55 20.98524
56 WB053 P -0.0174303 0.0055034 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 56 181.70564
57 WB053_3 G -1.3669331 0.0451303 Via & Conner 1995 24.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 57 22.15807
58 WB053_3 P 0.0035190 0.0048829 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 58 204.79462
59 WB053b G -1.3857924 0.0135062 Via & Conner 1995 80.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 59 74.04016
60 WB053b P -0.0149641 0.0060127 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 60 166.31377
61 WB053b_3 G -1.3855821 0.0116041 Via & Conner 1995 80.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 61 86.17632
62 WB053b_3 P 0.0045548 0.0064249 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 62 155.64431
63 WB057 G 1.1158349 0.0766307 Grill et al. 1997 24.0 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 63 13.04959
64 WB057 P 0.4308134 0.0117949 Grill et al. 1997 277.0 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 64 84.78246

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity estimates. Basically, these allow us to partition the variance into residual and between-study variance. We can also can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

    # Meta-analysis - Intercept Only Model
        # G 
            G_SMD_data <- subset(SMD_data, type == "G")
            
            model_intcp_SMD_het_G <- rma.mv(ES ~ 1, V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)

        # P 
            P_SMD_data <- subset(SMD_data, type == "P")
            
            model_intcp_SMD_het_P <- rma.mv(ES ~ 1, V = P_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = P_SMD_data)

    # Extract the variance estimates from the model
        #G
             stdy_sig_SMD_G <- model_intcp_SMD_het_G$sigma2[1]
              err_sig_SMD_G <- model_intcp_SMD_het_G$sigma2[2]

        #P
             stdy_sig_SMD_P <- model_intcp_SMD_het_P$sigma2[1]
              err_sig_SMD_P <- model_intcp_SMD_het_P$sigma2[2]

    # Calculate total sampling variance
        #G
                wi_SMD_G <- 1 / G_SMD_data$S_var
                vi_SMD_G <- sum(wi_SMD_G) * ( nrow(G_SMD_data) - 1) / ( (sum(wi_SMD_G)^2) - sum(wi_SMD_G))
                
        #P
                wi_SMD_P <- 1 / P_SMD_data$S_var
                vi_SMD_P <- sum(wi_SMD_P) * ( nrow(P_SMD_data) - 1) / ( (sum(wi_SMD_P)^2) - sum(wi_SMD_P))

    # Generate points estimates for I2 – heterogeneity at different levels
        #G
            I2total_SMD_G <- round((stdy_sig_SMD_G + err_sig_SMD_G) / (stdy_sig_SMD_G + err_sig_SMD_G + vi_SMD_G), digits = 2)

            I2stdy_SMD_G <- round((stdy_sig_SMD_G) / (stdy_sig_SMD_G + err_sig_SMD_G + vi_SMD_G), digits = 2)

        #P
            I2total_SMD_P <- round((stdy_sig_SMD_P + err_sig_SMD_P) / (stdy_sig_SMD_P + err_sig_SMD_P + vi_SMD_P), digits = 2)

            I2stdy_SMD_P <- round((stdy_sig_SMD_P) / (stdy_sig_SMD_P + err_sig_SMD_P + vi_SMD_P), digits = 2)

This will provide estimates of heterogeneity for G (I2total_SMD = 0.98 and between study heterogeneity as I2stdy_SMD = 0.84) and for P (I2total_SMD = 0.99 and between study heterogeneity as I2stdy_SMD = 0.97). While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should use a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the type of matrix (G or P), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actual effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this will reflect the median of the distribution of untransformed values rather than the mean.

    # Meta-regression models
        model_SMD_G <- metafor::rma.mv(ES ~ design + env.comp + scale(trait_num), V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)
        model_SMD_G_stress <- metafor::rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)
        model_SMD_G_half <- metafor::rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = G_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = G_SMD_data)
        
         n_design_SMD <- table(G_SMD_data$design)
        n_encComp_smd <- table(G_SMD_data$env.comp)
        
        # Marginalised means for different groups
        SMD_D_AN <- marginal(model_SMD_G, "design", G_SMD_data)
        SMD_D_BS <- marginal(model_SMD_G_stress, "design", G_SMD_data)

        AN_BS_SMD <- cbind(rbind(SMD_D_AN, SMD_D_BS), n_encComp_smd)

        SMD_D_full <- marginal(model_SMD_G, "env.comp", G_SMD_data)
        SMD_D_half <- marginal(model_SMD_G_half, "env.comp", G_SMD_data)

        design_SMD <- cbind(rbind(SMD_D_full, SMD_D_half), n_design_SMD)

        SMD_marg_table <- rbind(AN_BS_SMD, design_SMD)

        # We will just fit the P model
        model_SMD_P <- metafor::rma.mv(ES ~ design + env.comp + scale(trait_num), V = P_SMD_data$S_var, random = list(~1|stdy, ~1|err), data = P_SMD_data)

Once we have these models run we can have a look at the table of estimates:

    table_SMD <- data.frame(Est. = round(rbind(model_SMD_G$b, model_SMD_P$b), digits =3), L95CI = round(c(model_SMD_G$ci.lb, model_SMD_P$ci.lb), digits=3), U95CI = round(c(model_SMD_G$ci.ub, model_SMD_P$ci.ub), digits =3))
    rownames(table_SMD) <- c(paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_G"), paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_P"))
    kable(table_SMD) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept_G 0.623 0.086 1.160
Half-sib Breeding Design_G -1.092 -1.751 -0.434
Stressful Environment_G -0.011 -0.633 0.612
Trait Number_G -0.090 -0.359 0.179
Intercept_P 0.444 0.011 0.877
Half-sib Breeding Design_P -0.531 -1.059 -0.002
Stressful Environment_P -0.310 -0.809 0.190
Trait Number_P -0.200 -0.411 0.011

3.3 – Meta-analysis of Proportion of Variance Along max

First, we can have a look at the dataset more carefully before running the analysis.

    kable(var_along_logit_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")
X study type ES S_var authors year avg_n_type trait_num genus species design group env.comp species_full stdy err wi
1 RR073 G 0.3854686 0.0881714 Evans et al. 2015 66.0 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 1 11.341545
2 RR073 P 0.4469213 0.0303974 Evans et al. 2015 153.5 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 2 32.897531
3 RR076 G -0.9831790 0.1172049 Guan et al. 2016 69.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 3 8.532065
4 RR076 P -0.3773049 0.0031580 Guan et al. 2016 2525.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 4 316.654156
5 RR087a G 0.9952143 0.0739473 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 5 13.523150
6 RR087a P 0.2297012 0.0037580 Sae-Lim et al. 2013 2127.5 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 6 266.100828
7 RR087b G 0.4581764 0.0803353 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 7 12.447823
8 RR087b P -0.0555089 0.0031329 Sae-Lim et al. 2013 2577.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 8 319.193821
9 RR087c G -0.0984979 0.0750109 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 9 13.331403
10 RR087c P -0.5183770 0.0038306 Sae-Lim et al. 2013 2091.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 10 261.055378
11 sR097 G -1.1636097 0.3387223 Kasule 1991 26.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 11 2.952271
12 sR097 P -0.0071598 0.0253223 Kasule 1991 312.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 12 39.490913
13 sR110b G -2.8170762 0.2759173 Blanckenhorn & Heyland 2004 32.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 13 3.624274
14 sR110b P -0.5508976 0.0727283 Blanckenhorn & Heyland 2004 113.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 14 13.749807
15 sR110c G -1.0841059 0.2345270 Blanckenhorn & Heyland 2004 36.5 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 15 4.263901
16 sR110c P -0.6293073 0.0537564 Blanckenhorn & Heyland 2004 150.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 16 18.602440
17 sR115 G 0.3037840 0.0127205 Conner et al. 2003 471.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 17 78.613406
18 sR115 P 1.0542591 0.0051981 Conner et al. 2003 1010.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 18 192.378125
19 sR116 G -0.3731776 0.0201041 Delcourt et al. 2010 365.5 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 19 49.741126
20 sR116 P -0.0689952 0.0016041 Delcourt et al. 2010 3299.0 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 20 623.388723
21 sR145 G 0.5545398 0.0964286 Tucic et al. 1991 47.5 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 21 10.370367
22 sR145 P -0.1463109 0.0374098 Tucic et al. 1991 137.0 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 22 26.730963
23 WB002a G 0.0857084 0.1747969 Begin & Roff 2001 39.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 23 5.720927
24 WB002a P -0.2802787 0.0100423 Begin & Roff 2001 605.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 24 99.578399
25 WB003 G 0.5040949 0.0973839 Begin et al. 2004 59.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 25 10.268643
26 WB003 P 0.4346733 0.0086351 Begin et al. 2004 697.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 26 115.806945
27 WB003_2 G 0.0018268 0.0991137 Begin et al. 2004 58.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 27 10.089421
28 WB003_2 P 0.2619965 0.0090460 Begin et al. 2004 624.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 28 110.545546
29 WB003a G 0.2298057 0.1362611 Begin et al. 2004 43.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 29 7.338851
30 WB003a P 0.6737641 0.0191354 Begin et al. 2004 288.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 30 52.259231
31 WB003a_2 G 1.0096968 0.1659049 Begin et al. 2004 38.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 31 6.027551
32 WB003a_2 P 0.8927405 0.0305538 Begin et al. 2004 197.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 32 32.729101
33 WB007 G 0.7948070 0.0389207 Brock et al. 2010 131.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 33 25.693256
34 WB007 P 0.0222790 0.0089746 Brock et al. 2010 562.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 34 111.425492
35 WB010 G -0.1958279 0.0877272 Collins et al. 1999 101.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 35 11.398971
36 WB010 P -0.2305906 0.0103322 Collins et al. 1999 991.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 36 96.784861
37 WB011 G 1.4972372 0.0200357 Czesak & Fox 2003 404.0 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 37 49.910840
38 WB011 P 0.2038146 0.0029356 Czesak & Fox 2003 2667.0 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 38 340.646298
39 WB028 G 2.4763329 0.1339905 Klause & Morin 2001 60.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 39 7.463214
40 WB028 P 2.4941067 0.0130883 Klause & Morin 2001 600.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 40 76.404187
41 WB029 G 0.1310828 0.1289894 King et al. 2011 63.0 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 41 7.752577
42 WB029 P 0.1849028 0.0041783 King et al. 2011 1695.0 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 42 239.332622
43 WB033 G -0.1928236 0.1202423 Lau et al. 2014 60.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 43 8.316539
44 WB033 P 0.5169502 0.0128358 Lau et al. 2014 513.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 44 77.907073
45 WB033a G 0.7918906 0.1318105 Lau et al. 2014 60.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 45 7.586650
46 WB033a P 0.0045761 0.0141276 Lau et al. 2014 464.5 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 46 70.783494
47 WB035 G -1.4463889 0.0149326 Messina et al. 2013 432.0 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 47 66.967617
48 WB035 P -2.5735619 0.0030416 Messina et al. 2013 2204.0 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 48 328.770705
49 WB038 G 0.2309286 0.0439186 Rauter&Moore 2002 135.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 49 22.769383
50 WB038 P -0.5752013 0.0091143 Rauter&Moore 2002 639.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 50 109.717690
51 WB045 G 2.8750053 0.1531809 Simons & Roff 1996 39.0 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 51 6.528228
52 WB045 P 0.5532212 0.0149558 Simons & Roff 1996 371.0 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 52 66.863602
53 WB045a G 0.5798272 0.1588199 Simons & Roff 1996 39.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 53 6.296442
54 WB045a P 0.2880374 0.0152644 Simons & Roff 1996 381.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 54 65.511946
55 WB053 G -2.7022412 0.3643009 Via & Conner 1995 24.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 55 2.744983
56 WB053 P 0.2136209 0.0165484 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 56 60.428861
57 WB053_3 G -4.7140006 0.3802670 Via & Conner 1995 24.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 57 2.629731
58 WB053_3 P -0.2077135 0.0165950 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 58 60.259201
59 WB053b G -2.8298645 0.1024369 Via & Conner 1995 80.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 59 9.762108
60 WB053b P -1.0942523 0.0163428 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 60 61.189195
61 WB053b_3 G -4.7732524 0.1007585 Via & Conner 1995 80.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 61 9.924717
62 WB053b_3 P -0.1414758 0.0169504 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 62 58.995497
63 WB057 G -0.4329326 0.3313705 Grill et al. 1997 24.0 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 63 3.017770
64 WB057 P -0.5492158 0.0249628 Grill et al. 1997 277.0 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 64 40.059601

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity of effects. Basically, these allow us to partition the variance into residual, between-study variance and we can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

    # Multi-level meta-analysis
    # Subset the data and run the intercept only models
        # G
            G_varAlong_data <- subset(var_along_logit_data, type == "G")
            
            model_intcp_varAlong_het_G <- rma.mv(ES ~ 1, V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)

        # P
            P_varAlong_data <- subset(var_along_logit_data, type == "P")
            
            model_intcp_varAlong_het_P <- rma.mv(ES ~ 1, V = P_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = P_varAlong_data)

    # Extract the variance estimates from the model
        #G
            stdy_sig_varAlong_G <- model_intcp_varAlong_het_G$sigma2[1]
             err_sig_varAlong_G <- model_intcp_varAlong_het_G$sigma2[2]

        #P
             stdy_sig_varAlong_P <- model_intcp_varAlong_het_P$sigma2[1]
             err_sig_varAlong_P <- model_intcp_varAlong_het_P$sigma2[2]

    # Calculate total sampling variance
        #G
                wi_varAlong_G <- 1 / G_varAlong_data$S_var
                vi_varAlong_G <- sum(wi_varAlong_G) * ( nrow(G_varAlong_data) - 1) / ( (sum(wi_varAlong_G)^2) - sum(wi_varAlong_G))
        #P
                wi_varAlong_P <- 1 / P_varAlong_data$S_var
                vi_varAlong_P <- sum(wi_varAlong_P) * ( nrow(P_varAlong_data) - 1) / ( (sum(wi_varAlong_P)^2) - sum(wi_varAlong_P))

    # Generate points estimates for I2 – heterogenity at different levels
        #G
            I2total_varAlong_G <- round((stdy_sig_varAlong_G + err_sig_varAlong_G) / (stdy_sig_varAlong_G + err_sig_varAlong_G + vi_varAlong_G), digits = 2)

            I2stdy_varAlong_G <- round((stdy_sig_varAlong_G) / (stdy_sig_varAlong_G + err_sig_varAlong_G + vi_SMD_G), digits = 2)

        # P
            I2total_varAlong_P <- round((stdy_sig_varAlong_P + err_sig_varAlong_P) / (stdy_sig_varAlong_P + err_sig_varAlong_P + vi_varAlong_P), digits = 2)

            I2stdy_varAlong_P <- round((stdy_sig_varAlong_P) / (stdy_sig_varAlong_P + err_sig_varAlong_P + vi_SMD_P), digits = 2)

This will provide estimates of heterogeneity for G (I2total_Varlong = 0.97 and between study heterogeneity as I2stdy_varAlong = 0.73) and for P (I2total_varAlong = 0.99 and between study heterogeneity as I2stdy_varAlong = 0.82). While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the type of matrix (G or P), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actually effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this may reflect the median of the distribution of untransformed values rather than the mean. Models with a random slope for type within study converge nicely and so we include this in our models to account for correlations between G and P, not the least of which is the fact that they are based on the same underlying data.

   # Total variance change along g/p max.
        model_VarAlong_G <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)
        model_VarAlong_G_stress <- rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)
        model_VarAlong_G_half <- rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = G_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = G_varAlong_data)

        n_AN_S <- table(G_varAlong_data$env.comp)
        n_design <- table(G_varAlong_data$design)

        # Marginal means
        varalong_G_AN <- marginal(model_VarAlong_G, "design", data = G_varAlong_data)
        varalong_G_S  <- marginal(model_VarAlong_G_stress, "design", data = G_varAlong_data)

        env_cmp_varAlong <- rbind(varalong_G_AN, varalong_G_S)
        env_cmp_varAlong <- cbind(env_cmp_varAlong, n_AN_S)

        varalong_G_full <- marginal(model_VarAlong_G, "env.comp", G_varAlong_data)
        varalong_G_half <- marginal(model_VarAlong_G_half, "env.comp", G_varAlong_data)
          
        design_varAlong <- rbind(varalong_G_full, varalong_G_half)
        design_varAlong <- cbind(design_varAlong, n_design)

        varalong_marg_means <- rbind(env_cmp_varAlong, design_varAlong)

        model_VarAlong_P <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = P_varAlong_data$S_var, random = list(~1|stdy, ~1|err), data = P_varAlong_data)

Once we have these models run we can have a look at the table of estimates:

    table_varAlong <- data.frame(Est. = round(rbind(model_VarAlong_G$b, model_VarAlong_P$b), digits =3), L95CI = round(c(model_VarAlong_G$ci.lb, model_VarAlong_P$ci.lb), digits=3), U95CI = round(c(model_VarAlong_G$ci.ub, model_VarAlong_P$ci.ub), digits =3))
    rownames(table_varAlong) <- c(paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_G"), paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_P"))
    kable(table_varAlong) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept_G 0.466 -0.501 1.432
Half-sib Breeding Design_G -1.479 -2.681 -0.278
Stressful Environment_G 0.753 -0.387 1.893
Trait Number_G 0.453 -0.045 0.952
Intercept_P 0.320 -0.392 1.031
Half-sib Breeding Design_P -0.517 -1.393 0.360
Stressful Environment_P 0.011 -0.819 0.841
Trait Number_P 0.070 -0.289 0.430

3.4 – Meta-analysis of Angle Differences Across Environments

First, we can have a look at the dataset more carefully before running the analysis.

    kable(angle_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")
X study type ES S_var authors year avg_n_type trait_num genus species design group env.comp species_full stdy err wi
1 RR073 G 0.2525296 0.6071997 Evans et al. 2015 66.0 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 1 1.6469047
2 RR073 P 1.9626238 0.4822250 Evans et al. 2015 153.5 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 2 2.0737207
3 RR076 G -2.1552714 0.0548510 Guan et al. 2016 69.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 3 18.2311949
4 RR076 P -3.0958212 0.0199065 Guan et al. 2016 2525.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 4 50.2348513
5 RR087a G -2.4982323 0.3015269 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 5 3.3164541
6 RR087a P -1.8113762 0.0168979 Sae-Lim et al. 2013 2127.5 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 6 59.1788134
7 RR087b G -2.4223537 0.7267958 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 7 1.3759022
8 RR087b P -2.2142142 0.0399846 Sae-Lim et al. 2013 2577.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 8 25.0096430
9 RR087c G -2.4953079 0.7500006 Sae-Lim et al. 2013 100.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 9 1.3333322
10 RR087c P -1.4461472 0.0402833 Sae-Lim et al. 2013 2091.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 10 24.8242095
11 sR097 G -1.9352107 0.8905764 Kasule 1991 26.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 11 1.1228683
12 sR097 P -2.7588449 0.6487915 Kasule 1991 312.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 12 1.5413271
13 sR110b G -1.3485285 0.4753950 Blanckenhorn & Heyland 2004 32.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 13 2.1035140
14 sR110b P -1.5643303 0.1360035 Blanckenhorn & Heyland 2004 113.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 14 7.3527511
15 sR110c G -2.6295399 1.2333280 Blanckenhorn & Heyland 2004 36.5 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 15 0.8108143
16 sR110c P -2.2186783 0.2409692 Blanckenhorn & Heyland 2004 150.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 16 4.1499077
17 sR115 G -0.8877191 0.0241688 Conner et al. 2003 471.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 17 41.3756926
18 sR115 P -1.2472474 0.0135478 Conner et al. 2003 1010.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 18 73.8126450
19 sR116 G -3.4919763 0.2700806 Delcourt et al. 2010 365.5 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 19 3.7025989
20 sR116 P -4.8715979 0.1116487 Delcourt et al. 2010 3299.0 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 20 8.9566685
21 sR145 G 1.5171679 1.2131713 Tucic et al. 1991 47.5 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 21 0.8242859
22 sR145 P -2.2373713 0.1405605 Tucic et al. 1991 137.0 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 22 7.1143720
23 WB002a G -2.5297233 0.4065465 Begin & Roff 2001 39.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 23 2.4597430
24 WB002a P -2.9790205 0.1099553 Begin & Roff 2001 605.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 24 9.0946022
25 WB003 G -2.6296686 0.1720626 Begin et al. 2004 59.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 25 5.8118366
26 WB003 P -3.3843734 0.0551579 Begin et al. 2004 697.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 26 18.1297820
27 WB003_2 G -2.5908291 0.1617689 Begin et al. 2004 58.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 27 6.1816566
28 WB003_2 P -3.3986432 0.1291171 Begin et al. 2004 624.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 28 7.7449093
29 WB003a G -1.8234313 0.0993599 Begin et al. 2004 43.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 29 10.0644176
30 WB003a P -3.1296428 0.1030189 Begin et al. 2004 288.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 30 9.7069562
31 WB003a_2 G -2.1081164 0.0584278 Begin et al. 2004 38.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 31 17.1151325
32 WB003a_2 P -2.9168284 0.0972165 Begin et al. 2004 197.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 32 10.2863196
33 WB007 G 0.8783017 0.0321940 Brock et al. 2010 131.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 33 31.0616814
34 WB007 P -1.4142582 0.0086498 Brock et al. 2010 562.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 34 115.6091924
35 WB010 G -1.5891018 0.0166659 Collins et al. 1999 101.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 35 60.0028714
36 WB010 P -3.5814923 0.3460956 Collins et al. 1999 991.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 36 2.8893750
37 WB011 G -5.3388824 0.2657256 Czesak & Fox 2003 404.0 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 37 3.7632799
38 WB011 P -2.2309386 0.0172158 Czesak & Fox 2003 2667.0 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 38 58.0862199
39 WB028 G -0.2645351 0.0367901 Klause & Morin 2001 60.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 39 27.1812343
40 WB028 P -0.6601235 0.0217321 Klause & Morin 2001 600.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 40 46.0148668
41 WB029 G -4.3993865 0.6267768 King et al. 2011 63.0 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 41 1.5954643
42 WB029 P -0.8048412 0.0024736 King et al. 2011 1695.0 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 42 404.2631712
43 WB033 G -5.5539566 0.5108537 Lau et al. 2014 60.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 43 1.9575077
44 WB033 P -6.4522773 0.3105651 Lau et al. 2014 513.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 44 3.2199364
45 WB033a G -4.6843816 0.0960867 Lau et al. 2014 60.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 45 10.4072693
46 WB033a P -6.1899906 0.1131624 Lau et al. 2014 464.5 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 46 8.8368550
47 WB035 G -0.6864944 0.0175390 Messina et al. 2013 432.0 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 47 57.0157640
48 WB035 P -1.7223657 0.0215781 Messina et al. 2013 2204.0 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 48 46.3432925
49 WB038 G -2.0606936 0.2091481 Rauter&Moore 2002 135.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 49 4.7813013
50 WB038 P -2.7114595 0.0360297 Rauter&Moore 2002 639.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 50 27.7548998
51 WB045 G -0.2207482 0.0074963 Simons & Roff 1996 39.0 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 51 133.3985422
52 WB045 P -3.0747363 0.1228861 Simons & Roff 1996 371.0 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 52 8.1376165
53 WB045a G -2.8350612 0.2934883 Simons & Roff 1996 39.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 53 3.4072906
54 WB045a P -2.9119533 0.1649474 Simons & Roff 1996 381.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 54 6.0625388
55 WB053 G -2.8334129 0.4310894 Via & Conner 1995 24.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 55 2.3197045
56 WB053 P -3.5804857 1.1341314 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 56 0.8817320
57 WB053_3 G -2.5672827 1.1912009 Via & Conner 1995 24.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 57 0.8394889
58 WB053_3 P -0.4454588 0.0396640 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 58 25.2117556
59 WB053b G -3.0341318 0.1633039 Via & Conner 1995 80.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 59 6.1235516
60 WB053b P -1.6949102 0.0586645 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 60 17.0460785
61 WB053b_3 G -3.1392771 1.1214314 Via & Conner 1995 80.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 61 0.8917175
62 WB053b_3 P -3.0306846 0.2700453 Via & Conner 1995 480.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 62 3.7030825
63 WB057 G -1.9043333 0.1313720 Grill et al. 1997 24.0 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 63 7.6119738
64 WB057 P -3.1128513 0.0168762 Grill et al. 1997 277.0 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 64 59.2549553

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity estimates. Basically, these allow us to partition the variance into residual, between-study variance and we can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

    # Subset the data and run the intercept only models
     # G
        G_angle_data <- subset(angle_data, type == "G")
            
        model_intcp_angle_het_G <- rma.mv(ES ~ 1, V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)

    # P
        P_angle_data <- subset(angle_data, type == "P")
            
        model_intcp_angle_het_P <- rma.mv(ES ~ 1, V = P_angle_data$S_var, random = list(~1|stdy, ~1|err), data = P_angle_data)
            
    # Extract the variance estimates from the model     
        # G 
            stdy_sig_angle_G <- model_intcp_angle_het_G$sigma2[1]
             err_sig_angle_G <- model_intcp_angle_het_G$sigma2[2]
        # P
             stdy_sig_angle_P <- model_intcp_angle_het_P$sigma2[1]
             err_sig_angle_P <- model_intcp_angle_het_P$sigma2[2]

    # Calculate total sampling variance
        # G
                wi_angle_G <- 1 / G_angle_data$S_var
                vi_angle_G <- sum(wi_angle_G) * ( nrow(G_angle_data) - 1) / ( (sum(wi_angle_G)^2) - sum(wi_angle_G))
        # P 
                wi_angle_P <- 1 / P_angle_data$S_var
                vi_angle_P <- sum(wi_angle_P) * ( nrow(P_angle_data) - 1) / ( (sum(wi_angle_P)^2) - sum(wi_angle_P))
    
    # Generate points estimates for I2 – heterogenity at different levels
        #G
            I2total_angle_G <- round((stdy_sig_angle_G + err_sig_angle_G) / (stdy_sig_angle_G + err_sig_angle_G + vi_angle_G), digits = 2)

            I2stdy_angle_G <- round((stdy_sig_angle_G) / (stdy_sig_angle_G + err_sig_angle_G + vi_SMD_G), digits = 2)

        #P
            I2total_angle_P <- round((stdy_sig_angle_P + err_sig_angle_P) / (stdy_sig_angle_P + err_sig_angle_P + vi_angle_P), digits = 2)

            I2stdy_angle_P <- round((stdy_sig_angle_P) / (stdy_sig_angle_P + err_sig_angle_P + vi_SMD_P), digits = 2)

This will provide estimates of heterogeneity for G (I2total_angle = 0.98 and between study heterogeneity as I2stdy_angle = 0.85) and for P (I2total_angle = 0.99 and between study heterogeneity as I2stdy_angle = 0.88). While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the type of matrix (G or P), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actually effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this may reflect the median of the distribution of untransformed values rather than the mean. Models with a random slope for type within study converge nicely and so we include this in our models to account for correlations between G and P, not the least of which is the fact that they are based on the same underlying data.

    # Multi-level meta-regression models of the angle changes for G and P across environments. 
               model_Angle_G <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)
        model_Angle_G_stress <- rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)
          model_Angle_G_half <- rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = G_angle_data$S_var, random = list(~1|stdy, ~1|err), data = G_angle_data)

        n_AN_S_angle <- table(G_angle_data$env.comp)
        n_design_angle <- table(G_angle_data$design)

        # Marginal means
        angle_G_AN <- marginal(model_Angle_G, "design", data = G_angle_data)
        angle_G_S  <- marginal(model_Angle_G_stress, "design", data = G_angle_data)

        env_cmp_angle <- rbind(angle_G_AN, angle_G_S)
        env_cmp_angle <- cbind(env_cmp_angle, n_AN_S_angle)

        Angle_G_full <- marginal(model_Angle_G, "env.comp", G_angle_data)
        Angle_G_half <- marginal(model_Angle_G_half, "env.comp", G_angle_data)
          
        design_angle <- rbind(Angle_G_full, Angle_G_half)
        design_angle <- cbind(design_angle, n_design_angle)

        angle_marg_means <- rbind(env_cmp_angle, design_angle)

          model_Angle_P <- rma.mv(ES ~ design + env.comp + scale(trait_num), V = P_angle_data$S_var, random = list(~1|stdy, ~1|err), data = P_angle_data)

Once we have these models run we can have a look at the table of estimates:

    table_angle <- data.frame(Est. = round(rbind(model_Angle_G$b, model_Angle_P$b), digits =3), L95CI = round(c(model_Angle_G$ci.lb, model_Angle_P$ci.lb), digits=3), U95CI = round(c(model_Angle_G$ci.ub, model_Angle_P$ci.ub), digits =3))
    rownames(table_angle) <- c(paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_G"), paste0(c("Intercept", "Half-sib Breeding Design", "Stressful Environment", "Trait Number"), "_P"))
    kable(table_angle) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept_G -1.972 -3.015 -0.928
Half-sib Breeding Design_G -1.246 -2.551 0.058
Stressful Environment_G 1.513 0.273 2.753
Trait Number_G 0.891 0.353 1.428
Intercept_P -2.923 -4.062 -1.784
Half-sib Breeding Design_P -0.511 -1.920 0.898
Stressful Environment_P 1.711 0.369 3.052
Trait Number_P 0.097 -0.484 0.677

3.5 – Meta-analysis of Angle Differences Within Environments

First, we can have a look at the dataset more carefully before running the analysis.

    kable(angle_data_within) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")
X study env ES S_var environment authors year avg_n_env trait_num genus species design group env.comp species_full stdy err wi
1 RR073 B -0.3419286 0.5873615 E1 Evans et al. 2015 114.5 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 1 1.7025290
2 RR073 S -0.3158065 0.0463421 E2 Evans et al. 2015 105.0 8 Poecilia reticulata half-sib fish BS Poecilia_reticulata RR073 2 21.5786463
3 RR076 B -2.9218833 0.0317398 E1 Guan et al. 2016 1097.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 3 31.5061726
4 RR076 S -1.9911908 0.0618567 E2 Guan et al. 2016 1497.0 3 Scophthalmus maximus half-sib fish BS Scophthalmus_maximus RR076 4 16.1663968
5 RR087a A -0.9597197 0.1376774 E1 Sae-Lim et al. 2013 1232.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 5 7.2633537
6 RR087a N -1.4588856 0.0408052 E2 Sae-Lim et al. 2013 995.5 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 6 24.5066969
7 RR087b A -0.9584868 0.1239307 E1 Sae-Lim et al. 2013 1232.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 7 8.0690287
8 RR087b N -0.8841868 0.0427569 E2 Sae-Lim et al. 2013 1445.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 8 23.3880215
9 RR087c A -0.9148208 0.1567327 E1 Sae-Lim et al. 2013 1232.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 9 6.3802888
10 RR087c N -2.0659812 0.8070163 E2 Sae-Lim et al. 2013 959.0 3 Oncorhynchus mykiss full-sib fish AN Oncorhynchus_mykiss RR087 10 1.2391324
11 sR097 B -1.6290017 0.1391402 E1 Kasule 1991 169.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 11 7.1869948
12 sR097 S -2.8359926 1.3165859 E2 Kasule 1991 169.0 2 Dysdercus fasciatus half-sib insect BS Dysdercus_fasciatus sR097 12 0.7595403
13 sR110b B -3.2403778 0.7802790 E1 Blanckenhorn & Heyland 2004 82.5 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 13 1.2815928
14 sR110b S -2.9282182 1.4052676 E2 Blanckenhorn & Heyland 2004 62.5 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 14 0.7116082
15 sR110c B -3.7627363 1.2093555 E1 Blanckenhorn & Heyland 2004 105.5 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 15 0.8268867
16 sR110c S -3.2047958 1.3381308 E2 Blanckenhorn & Heyland 2004 81.0 2 Scathophaga stercoraria half-sib insect BS Scathophaga_stercoraria sR110 16 0.7473111
17 sR115 A -2.4518216 0.2393191 E1 Conner et al. 2003 741.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 17 4.1785219
18 sR115 N -0.1966421 0.0056296 E2 Conner et al. 2003 740.5 6 Raphanus raphanistrum half-sib plant AN Raphanus_raphanistrum sR115 18 177.6331442
19 sR116 A -2.0465360 0.0292682 E1 Delcourt et al. 2010 1840.5 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 19 34.1667396
20 sR116 N -1.8720304 0.0413366 E2 Delcourt et al. 2010 1824.0 8 Drosophila serrata half-sib insect AN Drosophila_serrata sR116 20 24.1916409
21 sR145 B 0.0908591 0.7744759 E1 Tucic et al. 1991 101.0 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 21 1.2911958
22 sR145 S 0.7859399 0.5287318 E2 Tucic et al. 1991 83.5 8 Acanthoscelides obtectus full-sib insect BS Acanthoscelides_obtectus sR145 22 1.8913179
23 WB002a A -1.8890462 0.2480546 E1 Begin & Roff 2001 372.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 23 4.0313701
24 WB002a N -1.8143146 0.1651411 E2 Begin & Roff 2001 272.0 5 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB002 24 6.0554263
25 WB003 A -2.9213749 0.1765009 E1 Begin et al. 2004 294.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 25 5.6656937
26 WB003 N -3.0281529 0.2354223 E2 Begin et al. 2004 387.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 26 4.2476864
27 WB003_2 A -2.9298980 0.1733582 E1 Begin et al. 2004 294.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 27 5.7684033
28 WB003_2 N -2.8386106 0.2638060 E2 Begin et al. 2004 462.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 28 3.7906645
29 WB003a A -2.6253466 0.1477390 E1 Begin et al. 2004 153.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 29 6.7686951
30 WB003a N -2.5919519 0.1156100 E2 Begin et al. 2004 82.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 30 8.6497698
31 WB003a_2 A -2.6548252 0.1429089 E1 Begin et al. 2004 153.5 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 31 6.9974646
32 WB003a_2 N -3.0880015 0.1029547 E2 Begin et al. 2004 178.0 5 Gryllus firmus full-sib insect AN Gryllus_firmus WB003 32 9.7130125
33 WB007 B 1.2753543 0.0083553 E1 Brock et al. 2010 366.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 33 119.6844085
34 WB007 S 2.4047661 0.1121431 E2 Brock et al. 2010 327.5 11 Brassica rapa half-sib plant BS Brassica_rapa WB007 34 8.9171751
35 WB010 A 5.5349758 1.2283841 E1 Collins et al. 1999 840.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 35 0.8140776
36 WB010 N 3.2749170 0.3578075 E2 Collins et al. 1999 252.5 4 Achroia grisella half-sib insect AN Achroia_grisella WB010 36 2.7947988
37 WB011 B 2.7850909 0.0260469 E1 Czesak & Fox 2003 1535.5 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 37 38.3922801
38 WB011 S 1.6471936 0.0029787 E2 Czesak & Fox 2003 1535.5 2 Stator limbatus half-sib insect BS Stator_limbatus WB011 38 335.7189771
39 WB028 B -2.3077811 0.6592608 E1 Klause & Morin 2001 330.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 39 1.5168505
40 WB028 S -4.5762373 1.0977522 E2 Klause & Morin 2001 330.0 2 Priophorus pallipes full-sib insect BS Priophorus_pallipes WB028 40 0.9109524
41 WB029 B 0.7016899 0.0028168 E1 King et al. 2011 862.5 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 41 355.0121381
42 WB029 S -0.6077572 0.0018738 E2 King et al. 2011 895.5 3 Gryllus firmus half-sib insect BS Gryllus_firmus WB029 42 533.6691304
43 WB033 A -4.6308526 0.0834496 E1 Lau et al. 2014 293.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 43 11.9832766
44 WB033 N -5.0518975 0.2464121 E2 Lau et al. 2014 280.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 44 4.0582423
45 WB033a A -4.1964803 0.0220943 E1 Lau et al. 2014 259.5 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 45 45.2604759
46 WB033a N -5.2362594 0.1008700 E2 Lau et al. 2014 265.0 4 Arabidopsis thaliana half-sib plant AN Arabidopsis_thaliana WB033 46 9.9137479
47 WB035 B -2.8244414 0.0259261 E1 Messina et al. 2013 1316.5 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 47 38.5711642
48 WB035 S -1.0705922 0.0404241 E2 Messina et al. 2013 1319.5 3 Callosobruchus maculatus half-sib insect BS Callosobruchus_maculatus WB035 48 24.7377434
49 WB038 B -3.0948590 0.2109398 E1 Rauter&Moore 2002 387.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 49 4.7406899
50 WB038 S -2.4697815 0.1683493 E2 Rauter&Moore 2002 387.0 4 Nicrophorus pustulatus half-sib insect BS Nicrophorus _pustulatus WB038 50 5.9400312
51 WB045 A -2.4725428 0.1708528 E1 Simons & Roff 1996 209.5 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 51 5.8529902
52 WB045 N -0.2165210 0.0005315 E2 Simons & Roff 1996 200.5 7 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 52 1881.2980685
53 WB045a A -2.2064629 0.2477518 E1 Simons & Roff 1996 215.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 53 4.0362980
54 WB045a N -3.0483883 0.2677812 E2 Simons & Roff 1996 205.0 6 Gryllus pennsylvanicus full-sib insect AN Gryllus_pennsylvanicus WB045 54 3.7343921
55 WB053 A 1.1219771 0.0184589 E1 Via & Conner 1995 252.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 55 54.1744813
56 WB053 N 0.6552539 0.0227371 E2 Via & Conner 1995 252.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 56 43.9809635
57 WB053_3 A 0.8219987 0.0188514 E1 Via & Conner 1995 252.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 57 53.0465318
58 WB053_3 N 1.7625772 0.7594263 E2 Via & Conner 1995 252.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 58 1.3167835
59 WB053b A 0.6628335 0.0030724 E1 Via & Conner 1995 280.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 59 325.4826934
60 WB053b N 1.2106564 0.0342508 E2 Via & Conner 1995 280.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 60 29.1963872
61 WB053b_3 A 0.4629305 0.0030871 E1 Via & Conner 1995 280.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 61 323.9316315
62 WB053b_3 N 0.4594170 0.0387154 E2 Via & Conner 1995 280.0 2 Tribolium castaneum half-sib insect AN Tribolium_castaneum WB053 62 25.8295272
63 WB057 A -2.7696908 0.1094535 E1 Grill et al. 1997 149.5 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 63 9.1362990
64 WB057 N -2.4674966 0.3017626 E2 Grill et al. 1997 151.5 4 Harmonia axyridis full-sib insect AN Harmonia_axyridis WB057 64 3.3138633

Once we are happy with the data we can run a simple intercept only meta-analytic model to estimate point estimates for heterogeneity estimates. Basically, these allow us to partition the variance into residual, between-study variance and we can estimate total sampling variance. All these measures are quite useful as they can tell us what might be driving differences in effects in the datasets.

     model_intcp_angle_within_het <- rma.mv(ES ~ 1, V = angle_data$S_var,   random = list(~1|stdy, ~1|err), data = angle_data_within)

            stdy_sig_angle_within <- model_intcp_angle_within_het$sigma2[1]
             err_sig_angle_within <- model_intcp_angle_within_het$sigma2[2]

        # Calculate total sampling variance
            wi_angle_within <- 1 / angle_data_within$S_var
            vi_angle_within <- sum(wi_angle_within) * ( nrow(angle_data_within) - 1) / ( (sum(wi_angle_within)^2) - sum(wi_angle_within))

        # Heterogeneity
    I2total_angle_within <- round((stdy_sig_angle_within + err_sig_angle_within) / (stdy_sig_angle_within + err_sig_angle_within + vi_angle_within), digits = 2)

    I2stdy_angle_within <- round((stdy_sig_angle_within) / (stdy_sig_angle_within + err_sig_angle_within + vi_angle_within), digits = 2)

Again, this will provide estimates of heterogeneity as the I2total_angle_within = 1 and between study heterogeneity as I2stdy_angle_within = 0.9. While these are useful, the model it self isn’t particularly helpful as we know, a priori that a number of moderator variables are likely driving effects. Hence, we should a meta-regression model. We explore only main effects of moderators especially given we don’t have a lot of data. We will account for differences in trait number across studies, the environment type of matrix (non-novel or novel), and the environmental comparison for the effect. In other words, whether the effect size comes from a study where stressful environmental conditions were indicated. We will run a model without an intercept for simplicity just in order to get estimates that we can back convert to get a handle on what the actually effect in the original units would be. Note that we are just doing some simple back transformations, but because of Jensen’s Inequality if the distributions are skewed this may reflect the median of the distribution of untransformed values rather than the mean. Models with a random slope for type within study converge nicely and so we include this in our models to account for correlations between G and P, not the least of which is the fact that they are based on the same underlying data.

    # Angle within environment multi-level meta regression models. Here we want to account for the fact that we have two environments per study and so these are correlated. 
        model_Angle_withinenv <- rma.mv(ES ~ environment + design + env.comp + scale(trait_num), V = angle_data_within$S_var, random = list(~1 + environment|stdy, ~1|err), data = angle_data_within)
            model_Angle_withinenv_stress <- rma.mv(ES ~ design + relevel(env.comp, ref = "BS") + scale(trait_num), V = angle_data_within$S_var, random = list(~1|stdy, ~1|err), data = angle_data_within)
          model_Angle_withinenv_half <- rma.mv(ES ~ relevel(design, ref = "half-sib") + env.comp + scale(trait_num), V = angle_data_within$S_var, random = list(~1|stdy, ~1|err), data = angle_data_within)

        n_AN_S_angle_within <- table(angle_data_within$env.comp)
        n_design_angle_within <- table(angle_data_within$design)

        # Marginal means
        angle_within_G_AN <- marginal(model_Angle_withinenv, "design", data = angle_data_within)
        angle_within_G_S  <- marginal(model_Angle_withinenv_stress, "design", data = angle_data_within)

        env_cmp_angle_within <- rbind(angle_within_G_AN, angle_within_G_S)
        env_cmp_angle_within <- cbind(env_cmp_angle_within, n_AN_S_angle_within)

        Angle_within_G_full <- marginal(model_Angle_withinenv, "env.comp", angle_data_within)
        Angle_within_G_half <- marginal(model_Angle_withinenv_half, "env.comp", angle_data_within)
          
        design_angle_within <- rbind(Angle_within_G_full, Angle_within_G_half)
        design_angle_within <- cbind(design_angle_within, n_design_angle_within)

        angle_within_marg_means <- rbind(env_cmp_angle_within, design_angle_within)

Once we have these models run we can have a look at the table of estimates:

    table_angle_within <- data.frame(Est. = round(model_Angle_withinenv$b, digits =3), L95CI = round(model_Angle_withinenv$ci.lb, digits=3), U95CI = round(model_Angle_withinenv$ci.ub, digits =3))
    rownames(table_angle_within) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Stressful Environment", "Trait Number")

    kable(table_angle_within) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept -2.130 -3.858 -0.402
Novel Environment 0.147 -0.305 0.599
Half-sib Breeding Design 1.078 -1.024 3.180
Stressful Environment 0.095 -1.886 2.076
Trait Number 0.754 -0.050 1.558

4.0 – Meta-analysis 2: Evolvability of environmentally induced phenotypes

4.1 – Load Data

              evolve_max_data <- read.csv("./data_evolvability/evolve_max_data_standardised.csv")
         evolve_max_data_test <- read.csv("./data_evolvability/evolve_max_data_test.csv")
            evolve_angle_data <- read.csv("./data_evolvability/evolve_angle_data_standardised.csv")

            evolve_max_data$stress <- ifelse(evolve_max_data$env == "S", "S", "A")
            evolve_angle_data$stress <- ifelse(evolve_angle_data$env == "S", "S", "A")

4.2 – Meta-analysis of \(\pi_{e}\)

Again, we can view the raw data as normal.

    kable(evolve_max_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")
authors year study env type n trait_num genus species design group name.type name.env ES S_var err environment stdy
Beacham 1988 RR062 N G 10 3 Oncorhynchus gorbuscha half-sib fish RR062.G RR062.N 0.5428624 0.1240245 1 Novel RR062
Beacham 1988 RR062_b N G 10 3 Oncorhynchus gorbuscha half-sib fish RR062_b.G RR062_b.N 0.2101515 1.7883816 2 Novel RR062
Beacham 1988 RR062b N G 10 3 Oncorhynchus keta half-sib fish RR062b.G RR062b.N -1.2101635 0.8238148 3 Novel RR062
Beacham 1988 RR062b_b N G 10 3 Oncorhynchus keta half-sib fish RR062b_b.G RR062b_b.N -4.5066487 0.2344635 4 Novel RR062
Evans et al. 2015 RR073 S G 64 8 Poecilia reticulata half-sib fish RR073.G RR073.S 0.4505784 0.0172563 5 Novel RR073
Guan et al. 2016 RR076 S G 69 3 Scophthalmus maximus half-sib fish RR076.G RR076.S 1.3571034 0.0152012 6 Novel RR076
Sae-Lim et al. 2013 RR087a N G 100 3 Oncorhynchus mykiss full-sib fish RR087a.G RR087a.N 0.6551657 0.0294507 7 Novel RR087
Sae-Lim et al. 2013 RR087b N G 100 3 Oncorhynchus mykiss full-sib fish RR087b.G RR087b.N -0.1804580 0.0015965 8 Novel RR087
Sae-Lim et al. 2013 RR087c N G 100 3 Oncorhynchus mykiss full-sib fish RR087c.G RR087c.N 0.6335908 0.1395597 9 Novel RR087
Kasule 1991 sR097 S G 26 2 Dysdercus fasciatus half-sib insect sR097.G sR097.S 5.0286449 4.4480211 10 Novel sR097
Blanckenhorn & Heyland 2004 sR110b S G 27 2 Scathophaga stercoraria half-sib insect sR110b.G sR110b.S 1.7989274 0.6511523 11 Novel sR110
Blanckenhorn & Heyland 2004 sR110c S G 33 2 Scathophaga stercoraria half-sib insect sR110c.G sR110c.S 1.4345503 0.2583591 12 Novel sR110
Conner et al. 2003 sR115 N G 593 6 Raphanus raphanistrum half-sib plant sR115.G sR115.N 0.8059430 0.0037772 13 Novel sR115
Delcourt et al. 2010 sR116 N G 365 8 Drosophila serrata half-sib insect sR116.G sR116.N -4.0826705 0.0077509 14 Novel sR116
Tucic et al. 1991 sR145 S G 43 8 Acanthoscelides obtectus full-sib insect sR145.G sR145.S 0.0775910 0.0771394 15 Novel sR145
Auld 2010 WB001 S G 30 3 Physa acuta full-sib mollusc WB001.G WB001.S 0.2444737 0.6876215 16 Novel WB001
Auld 2010 WB001_b S G 30 3 Physa acuta full-sib mollusc WB001_b.G WB001_b.S 2.5420243 0.7899251 17 Novel WB001
Auld 2010 WB001_c S G 30 3 Physa acuta full-sib mollusc WB001_c.G WB001_c.S 1.3910750 0.1019580 18 Novel WB001
Begin & Roff 2001 WB002a N G 39 5 Gryllus pennsylvanicus full-sib insect WB002a.G WB002a.N 0.9910524 0.0537902 19 Novel WB002
Begin et al. 2004 WB003 N G 62 5 Gryllus firmus full-sib insect WB003.G WB003.N -3.0566626 0.0450678 20 Novel WB003
Begin et al. 2004 WB003_2 N G 60 5 Gryllus firmus full-sib insect WB003_2.G WB003_2.N -2.3639847 0.0373027 21 Novel WB003
Begin et al. 2004 WB003a N G 41 5 Gryllus firmus full-sib insect WB003a.G WB003a.N 0.8568436 0.0204842 22 Novel WB003
Begin et al. 2004 WB003a_2 N G 32 5 Gryllus firmus full-sib insect WB003a_2.G WB003a_2.N -0.4183020 0.0048221 23 Novel WB003
Brock et al. 2010 WB007 S G 131 11 Brassica rapa half-sib plant WB007.G WB007.S -3.1248238 0.0085661 24 Novel WB007
Collins et al. 1999 WB010 N G 62 4 Achroia grisella half-sib insect WB010.G WB010.N 1.1622610 0.0075125 25 Novel WB010
Czesak & Fox 2003 WB011 S G 404 2 Stator limbatus half-sib insect WB011.G WB011.S -2.9885446 0.0001467 26 Novel WB011
Delcourt et al. 2009 WB014 N G 317 2 Drosophila serrata half-sib insect WB014.G WB014.N 5.2775656 0.0003394 27 Novel WB014
Engqvist 2007 WB016 S G 27 2 Panorpa cognata full-sib insect WB016.G WB016.S 10.0864227 2.2858387 28 Novel WB016
Guntrip et al. 1997 WB022 N G 40 2 Callosobruchus maculatus half-sib insect WB022.G WB022.N 0.8558881 0.1851336 29 Novel WB022
Holloway et al. 1990 WB025 N G 69 4 Sitophilus oryzae half-sib insect WB025.G WB025.N 1.2765941 0.0699061 30 Novel WB025
Klause & Morin 2001 WB028 S G 60 2 Priophorus pallipes full-sib insect WB028.G WB028.S -1.1091105 0.0060447 31 Novel WB028
King et al. 2011 WB029 S G 63 3 Gryllus firmus half-sib insect WB029.G WB029.S -2.9537509 0.0131215 32 Novel WB029
Lau et al. 2014 WB033 N G 60 4 Arabidopsis thaliana half-sib plant WB033.G WB033.N 1.0437744 0.0000261 33 Novel WB033
Lau et al. 2014 WB033a N G 60 4 Arabidopsis thaliana half-sib plant WB033a.G WB033a.N -2.9124026 0.0000010 34 Novel WB033
Messina et al. 2013 WB035 S G 432 3 Callosobruchus maculatus half-sib insect WB035.G WB035.S 5.5948655 0.7656785 35 Novel WB035
Paccard et al. 2013 WB036 N G 22 7 Arabidopsis lyrata full-sib plant WB036.G WB036.N 0.7126878 0.0132725 36 Novel WB036
Paccard et al. 2013 WB036A N G 22 7 Arabidopsis lyrata full-sib plant WB036A.G WB036A.N 1.8722282 0.0597320 37 Novel WB036
Rauter&Moore 2002 WB038 S G 135 4 Nicrophorus pustulatus half-sib insect WB038.G WB038.S 3.4314013 0.0984880 38 Novel WB038
Relyea 2005 WB039 S G 21 9 Rana sylvatica half-sib amphibian WB039.G WB039.S -0.1582131 0.1651519 39 Novel WB039
Sherrard et al. 2009 WB044 S G 26 10 Avena barbata half-sib plant WB044.G WB044.S -3.8866656 0.1485572 40 Novel WB044
Simons & Roff 1996 WB045 N G 39 7 Gryllus pennsylvanicus full-sib insect WB045.G WB045.N 0.5875816 0.0024073 41 Novel WB045
Simons & Roff 1996 WB045a N G 39 6 Gryllus pennsylvanicus full-sib insect WB045a.G WB045a.N 3.2819619 0.0732472 42 Novel WB045
Simonsen & Stinchcombe 2001 WB046 S G 50 3 Ipomoea hederacea half-sib plant WB046.G WB046.S 4.1499697 0.9443279 43 Novel WB046
Via & Conner 1995 WB053 N G 24 2 Tribolium castaneum half-sib insect WB053.G WB053.N 0.6486913 0.0301427 44 Novel WB053
Via & Conner 1995 WB053_3 N G 24 2 Tribolium castaneum half-sib insect WB053_3.G WB053_3.N 0.0056850 0.2296783 45 Novel WB053
Via & Conner 1995 WB053b N G 80 2 Tribolium castaneum half-sib insect WB053b.G WB053b.N 0.7647140 0.0100357 46 Novel WB053
Via & Conner 1995 WB053b_3 N G 80 2 Tribolium castaneum half-sib insect WB053b_3.G WB053b_3.N -0.5880988 0.0682137 47 Novel WB053
Service & Rose 1985 WB056 N G 64 2 Drosophila melanogaster half-sib insect WB056.G WB056.N 1.8715821 0.4847043 48 Novel WB056
Grill et al. 1997 WB057 N G 26 4 Harmonia axyridis full-sib insect WB057.G WB057.N -1.3737456 0.0503796 49 Novel WB057
Beacham 1988 RR062 A G 10 3 Oncorhynchus gorbuscha half-sib fish RR062.G RR062.A -1.8588986 0.4660717 50 Non-novel RR062
Beacham 1988 RR062_b A G 10 3 Oncorhynchus gorbuscha half-sib fish RR062_b.G RR062_b.A -0.4636018 0.6964657 51 Non-novel RR062
Beacham 1988 RR062b A G 10 3 Oncorhynchus keta half-sib fish RR062b.G RR062b.A -2.4161392 0.4171573 52 Non-novel RR062
Beacham 1988 RR062b_b A G 10 3 Oncorhynchus keta half-sib fish RR062b_b.G RR062b_b.A -0.2500228 0.5239350 53 Non-novel RR062
Evans et al. 2015 RR073 B G 68 8 Poecilia reticulata half-sib fish RR073.G RR073.B 1.1675329 0.0670842 54 Non-novel RR073
Guan et al. 2016 RR076 B G 69 3 Scophthalmus maximus half-sib fish RR076.G RR076.B 1.0152374 0.0082485 55 Non-novel RR076
Sae-Lim et al. 2013 RR087a A G 100 3 Oncorhynchus mykiss full-sib fish RR087a.G RR087a.A 1.3459910 0.1654308 56 Non-novel RR087
Sae-Lim et al. 2013 RR087b A G 100 3 Oncorhynchus mykiss full-sib fish RR087b.G RR087b.A -0.0889492 0.0226032 57 Non-novel RR087
Sae-Lim et al. 2013 RR087c A G 100 3 Oncorhynchus mykiss full-sib fish RR087c.G RR087c.A 0.5866717 0.1139000 58 Non-novel RR087
Kasule 1991 sR097 B G 26 2 Dysdercus fasciatus half-sib insect sR097.G sR097.B 4.3263515 1.7230775 59 Non-novel sR097
Blanckenhorn & Heyland 2004 sR110b B G 37 2 Scathophaga stercoraria half-sib insect sR110b.G sR110b.B -0.1226381 0.0067473 60 Non-novel sR110
Blanckenhorn & Heyland 2004 sR110c B G 40 2 Scathophaga stercoraria half-sib insect sR110c.G sR110c.B 0.6388066 0.0367883 61 Non-novel sR110
Conner et al. 2003 sR115 A G 350 6 Raphanus raphanistrum half-sib plant sR115.G sR115.A 3.4973585 0.1625815 62 Non-novel sR115
Delcourt et al. 2010 sR116 A G 366 8 Drosophila serrata half-sib insect sR116.G sR116.A -4.4340787 0.0079699 63 Non-novel sR116
Tucic et al. 1991 sR145 B G 52 8 Acanthoscelides obtectus full-sib insect sR145.G sR145.B 0.3300179 0.1013896 64 Non-novel sR145
Auld 2010 WB001 B G 30 3 Physa acuta full-sib mollusc WB001.G WB001.B -1.2919352 0.2238658 65 Non-novel WB001
Auld 2010 WB001_b B G 30 3 Physa acuta full-sib mollusc WB001_b.G WB001_b.B 1.9766137 0.7477113 66 Non-novel WB001
Auld 2010 WB001_c B G 30 3 Physa acuta full-sib mollusc WB001_c.G WB001_c.B 1.1365886 0.6238066 67 Non-novel WB001
Begin & Roff 2001 WB002a A G 39 5 Gryllus pennsylvanicus full-sib insect WB002a.G WB002a.A 0.9551444 0.0507751 68 Non-novel WB002
Begin et al. 2004 WB003 A G 56 5 Gryllus firmus full-sib insect WB003.G WB003.A -3.0254976 0.0642904 69 Non-novel WB003
Begin et al. 2004 WB003_2 A G 56 5 Gryllus firmus full-sib insect WB003_2.G WB003_2.A -2.5292244 0.0540068 70 Non-novel WB003
Begin et al. 2004 WB003a A G 45 5 Gryllus firmus full-sib insect WB003a.G WB003a.A 0.4031230 0.0321589 71 Non-novel WB003
Begin et al. 2004 WB003a_2 A G 45 5 Gryllus firmus full-sib insect WB003a_2.G WB003a_2.A -0.5568401 0.0224074 72 Non-novel WB003
Brock et al. 2010 WB007 B G 132 11 Brassica rapa half-sib plant WB007.G WB007.B -0.9073448 0.0187114 73 Non-novel WB007
Collins et al. 1999 WB010 A G 141 4 Achroia grisella half-sib insect WB010.G WB010.A 1.4115091 0.0019616 74 Non-novel WB010
Czesak & Fox 2003 WB011 B G 404 2 Stator limbatus half-sib insect WB011.G WB011.B -3.0142257 0.0006719 75 Non-novel WB011
Delcourt et al. 2009 WB014 A G 319 2 Drosophila serrata half-sib insect WB014.G WB014.A 5.2246319 0.0017553 76 Non-novel WB014
Engqvist 2007 WB016 B G 27 2 Panorpa cognata full-sib insect WB016.G WB016.B 12.5280490 4.9242628 77 Non-novel WB016
Guntrip et al. 1997 WB022 A G 40 2 Callosobruchus maculatus half-sib insect WB022.G WB022.A -1.3775179 0.0789878 78 Non-novel WB022
Holloway et al. 1990 WB025 A G 65 4 Sitophilus oryzae half-sib insect WB025.G WB025.A 2.1329860 0.0104086 79 Non-novel WB025
Klause & Morin 2001 WB028 B G 60 2 Priophorus pallipes full-sib insect WB028.G WB028.B -1.2711403 0.0982403 80 Non-novel WB028
King et al. 2011 WB029 B G 63 3 Gryllus firmus half-sib insect WB029.G WB029.B -2.7665833 0.0132234 81 Non-novel WB029
Lau et al. 2014 WB033 A G 60 4 Arabidopsis thaliana half-sib plant WB033.G WB033.A 1.0449722 0.0000227 82 Non-novel WB033
Lau et al. 2014 WB033a A G 60 4 Arabidopsis thaliana half-sib plant WB033a.G WB033a.A -2.9069230 0.0000013 83 Non-novel WB033
Messina et al. 2013 WB035 B G 432 3 Callosobruchus maculatus half-sib insect WB035.G WB035.B 0.9121341 0.0020951 84 Non-novel WB035
Paccard et al. 2013 WB036 A G 22 7 Arabidopsis lyrata full-sib plant WB036.G WB036.A 0.7286765 0.1116486 85 Non-novel WB036
Paccard et al. 2013 WB036A A G 22 7 Arabidopsis lyrata full-sib plant WB036A.G WB036A.A 2.3949669 0.0263805 86 Non-novel WB036
Rauter&Moore 2002 WB038 B G 135 4 Nicrophorus pustulatus half-sib insect WB038.G WB038.B 3.7426764 0.1518899 87 Non-novel WB038
Relyea 2005 WB039 B G 21 9 Rana sylvatica half-sib amphibian WB039.G WB039.B -0.3485422 0.1641508 88 Non-novel WB039
Sherrard et al. 2009 WB044 B G 26 10 Avena barbata half-sib plant WB044.G WB044.B -4.1224059 0.1631106 89 Non-novel WB044
Simons & Roff 1996 WB045 A G 39 7 Gryllus pennsylvanicus full-sib insect WB045.G WB045.A 0.1642072 0.0270893 90 Non-novel WB045
Simons & Roff 1996 WB045a A G 39 6 Gryllus pennsylvanicus full-sib insect WB045a.G WB045a.A 3.2745521 0.2456119 91 Non-novel WB045
Simonsen & Stinchcombe 2001 WB046 B G 50 3 Ipomoea hederacea half-sib plant WB046.G WB046.B 2.0149730 0.2862508 92 Non-novel WB046
Via & Conner 1995 WB053 A G 24 2 Tribolium castaneum half-sib insect WB053.G WB053.A 1.0237109 0.0023335 93 Non-novel WB053
Via & Conner 1995 WB053_3 A G 24 2 Tribolium castaneum half-sib insect WB053_3.G WB053_3.A -1.1432159 0.0018889 94 Non-novel WB053
Via & Conner 1995 WB053b A G 80 2 Tribolium castaneum half-sib insect WB053b.G WB053b.A 1.0418965 0.0006471 95 Non-novel WB053
Via & Conner 1995 WB053b_3 A G 80 2 Tribolium castaneum half-sib insect WB053b_3.G WB053b_3.A -1.1265929 0.0005348 96 Non-novel WB053
Service & Rose 1985 WB056 A G 64 2 Drosophila melanogaster half-sib insect WB056.G WB056.A 3.4838459 0.0562978 97 Non-novel WB056
Grill et al. 1997 WB057 A G 22 4 Harmonia axyridis full-sib insect WB057.G WB057.A -0.6138681 0.0260655 98 Non-novel WB057

We can now estimate heterogeneity using a multi-level meta-analytic model

    mod_evolve_max_int <-  rma.mv(ES ~ 1, V = evolve_max_data$S_var, random = list(~1|stdy, ~1|err), data = evolve_max_data)

            stdy_evolve_max_int <- mod_evolve_max_int$sigma2[1]
             err_evolve_max_int <- mod_evolve_max_int$sigma2[2]
    
    # Calculate total sampling variance
            wi_evolve_max <- 1 / evolve_max_data$S_var
            vi_evolve_max <- sum(wi_evolve_max) * ( nrow(evolve_max_data) - 1) / ( (sum(wi_evolve_max)^2) - sum(wi_evolve_max))

            I2total_evolve_max_data <- round((stdy_evolve_max_int + err_evolve_max_int) / (stdy_evolve_max_int + err_evolve_max_int + vi_evolve_max), digits = 2)

            I2stdy_evolve_max_data <- round((stdy_evolve_max_int) / (stdy_evolve_max_int + err_evolve_max_int + vi_evolve_max), digits = 2)

OK, but again, we have a priori moderators we think should explain these

            mod_evolve_max_trait_env <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)
            mod_evolve_max_trait_env_stress <-  rma.mv(ES ~ 1 + environment + design + stress + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)

        # calculate marginal mean for non-novel, marginalised over design type. First, refit model
            mod_evolve_max_trait_env_novel <-  rma.mv(ES ~ 1 + relevel(environment, ref = "Novel") + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)          
            
            marinal_means_evolve_max_A_N <- rbind(marginal(mod_evolve_max_trait_env, "design", data = evolve_max_data), marginal(mod_evolve_max_trait_env_novel, "design", data = evolve_max_data))
            n_A_N <- table(evolve_max_data$environment)
            marinal_means_evolve_max_A_N <- cbind(marinal_means_evolve_max_A_N, n_A_N)
            
        # Now we can do this over design
            mod_evolve_max_trait_env_half <-  rma.mv(ES ~ 1 + environment + relevel(design, ref = "half-sib") + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)
            marinal_means_evolve_max_design <- rbind(marginal(mod_evolve_max_trait_env, "environment", data = evolve_max_data), marginal(mod_evolve_max_trait_env_half, "environment", data = evolve_max_data))
            n_design <- table(evolve_max_data$design)
            marinal_means_evolve_max_design <- cbind(marinal_means_evolve_max_design, n_design)

            marg_evolve_max <- rbind(marinal_means_evolve_max_A_N, marinal_means_evolve_max_design)

Once we have these models run we can have a look at the table of estimates:

    table_max_evolve <- data.frame(Est. = round(mod_evolve_max_trait_env$b, digits =3), L95CI = round(mod_evolve_max_trait_env$ci.lb, digits=3), U95CI = round(mod_evolve_max_trait_env$ci.ub, digits =3))
    rownames(table_max_evolve) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Trait Number")   
    kable(table_max_evolve) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept 1.268 -0.356 2.893
Novel Environment 0.092 -0.417 0.601
Half-sib Breeding Design -0.754 -2.685 1.176
Trait Number -1.073 -1.873 -0.274

To test whether there is more additive genetic variation along the plasticity vector than expected we also run the model for the difference between \(\pi_{e}\) and \(\pi_{0}\).

        mod_evolve_max_test_env <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_max_data_test$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data_test)

Once we have this model run we can have a look at the table of estimates:

4.3 – Table S1

Table S1 – Estimates testing whether the angle of gmax is aligned more than simply by chance. We are specifically interested in the intercept here and when the 95% confidence intervals exclude zero, there is evidence that alignment between plasticity vector and \(g_{max}\) is more then what we expect by chance alone.

    table_max_evolve_test <- data.frame(Est. = round(mod_evolve_max_test_env$b, digits =3), L95CI = round(mod_evolve_max_test_env$ci.lb, digits=3), U95CI = round(mod_evolve_max_test_env$ci.ub, digits =3))
    rownames(table_max_evolve_test) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Trait Number")  
    kable(table_max_evolve_test) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept 1.108 0.393 1.822
Novel Environment 0.040 -0.241 0.321
Half-sib Breeding Design -0.573 -1.420 0.274
Trait Number 0.164 -0.192 0.521

4.4 – Meta-analysis of \(\theta_{e}\)

In case a reader is interested, we’ll still spit out the raw data:

    kable(evolve_angle_data) %>%
    kable_styling(font = 10) %>%
    scroll_box(width = "100%", height = "400px")
authors year study env type n trait_num genus species design group name.type name.env ES S_var err environment stdy
Beacham 1988 RR062 N G 10 3 Oncorhynchus gorbuscha half-sib fish RR062.G RR062.N -0.2804839 0.0692202 1 Novel RR062
Beacham 1988 RR062_b N G 10 3 Oncorhynchus gorbuscha half-sib fish RR062_b.G RR062_b.N 1.0696574 2.1866890 2 Novel RR062
Beacham 1988 RR062b N G 10 3 Oncorhynchus keta half-sib fish RR062b.G RR062b.N 2.1656787 1.4562036 3 Novel RR062
Beacham 1988 RR062b_b N G 10 3 Oncorhynchus keta half-sib fish RR062b_b.G RR062b_b.N 3.0168584 0.2128005 4 Novel RR062
Evans et al. 2015 RR073 S G 64 8 Poecilia reticulata half-sib fish RR073.G RR073.S -0.2136747 0.0122647 5 Novel RR073
Guan et al. 2016 RR076 S G 69 3 Scophthalmus maximus half-sib fish RR076.G RR076.S -0.8207908 0.0063809 6 Novel RR076
Sae-Lim et al. 2013 RR087a N G 100 3 Oncorhynchus mykiss full-sib fish RR087a.G RR087a.N -0.2667054 0.0127446 7 Novel RR087
Sae-Lim et al. 2013 RR087b N G 100 3 Oncorhynchus mykiss full-sib fish RR087b.G RR087b.N 0.1299726 0.0008169 8 Novel RR087
Sae-Lim et al. 2013 RR087c N G 100 3 Oncorhynchus mykiss full-sib fish RR087c.G RR087c.N 0.2813025 0.1212850 9 Novel RR087
Kasule 1991 sR097 S G 26 2 Dysdercus fasciatus half-sib insect sR097.G sR097.S -2.7552460 1.2302908 10 Novel sR097
Blanckenhorn & Heyland 2004 sR110b S G 27 2 Scathophaga stercoraria half-sib insect sR110b.G sR110b.S -0.8894004 0.2351409 11 Novel sR110
Blanckenhorn & Heyland 2004 sR110c S G 33 2 Scathophaga stercoraria half-sib insect sR110c.G sR110c.S -0.7073991 0.1014419 12 Novel sR110
Conner et al. 2003 sR115 N G 593 6 Raphanus raphanistrum half-sib plant sR115.G sR115.N -0.4036144 0.0043551 13 Novel sR115
Delcourt et al. 2010 sR116 N G 365 8 Drosophila serrata half-sib insect sR116.G sR116.N 3.0267427 0.0073985 14 Novel sR116
Tucic et al. 1991 sR145 S G 43 8 Acanthoscelides obtectus full-sib insect sR145.G sR145.S 0.2679014 0.1688566 15 Novel sR145
Auld 2010 WB001 S G 30 3 Physa acuta full-sib mollusc WB001.G WB001.S 1.6941891 1.8150091 16 Novel WB001
Auld 2010 WB001_b S G 30 3 Physa acuta full-sib mollusc WB001_b.G WB001_b.S -1.4002198 0.2552470 17 Novel WB001
Auld 2010 WB001_c S G 30 3 Physa acuta full-sib mollusc WB001_c.G WB001_c.S -0.7901196 0.0400435 18 Novel WB001
Begin & Roff 2001 WB002a N G 39 5 Gryllus pennsylvanicus full-sib insect WB002a.G WB002a.N -0.5388656 0.0264206 19 Novel WB002
Begin et al. 2004 WB003 N G 62 5 Gryllus firmus full-sib insect WB003.G WB003.N 2.4367906 0.0410906 20 Novel WB003
Begin et al. 2004 WB003_2 N G 60 5 Gryllus firmus full-sib insect WB003_2.G WB003_2.N 1.8226200 0.0232361 21 Novel WB003
Begin et al. 2004 WB003a N G 41 5 Gryllus firmus full-sib insect WB003a.G WB003a.N -0.5049055 0.0085218 22 Novel WB003
Begin et al. 2004 WB003a_2 N G 32 5 Gryllus firmus full-sib insect WB003a_2.G WB003a_2.N 0.2796241 0.0021071 23 Novel WB003
Brock et al. 2010 WB007 S G 131 11 Brassica rapa half-sib plant WB007.G WB007.S 2.0510321 0.0038105 24 Novel WB007
Collins et al. 1999 WB010 N G 62 4 Achroia grisella half-sib insect WB010.G WB010.N -0.7145406 0.0036595 25 Novel WB010
Czesak & Fox 2003 WB011 S G 404 2 Stator limbatus half-sib insect WB011.G WB011.S 1.8190447 0.0000475 26 Novel WB011
Delcourt et al. 2009 WB014 N G 317 2 Drosophila serrata half-sib insect WB014.G WB014.N -3.0455236 0.0000925 27 Novel WB014
Engqvist 2007 WB016 S G 27 2 Panorpa cognata full-sib insect WB016.G WB016.S -5.4896469 0.5746361 28 Novel WB016
Guntrip et al. 1997 WB022 N G 40 2 Callosobruchus maculatus half-sib insect WB022.G WB022.N -0.2539235 0.0864514 29 Novel WB022
Holloway et al. 1990 WB025 N G 69 4 Sitophilus oryzae half-sib insect WB025.G WB025.N -0.6432784 0.0496260 30 Novel WB025
Klause & Morin 2001 WB028 S G 60 2 Priophorus pallipes full-sib insect WB028.G WB028.S 0.7417540 0.0024437 31 Novel WB028
King et al. 2011 WB029 S G 63 3 Gryllus firmus half-sib insect WB029.G WB029.S 1.9110061 0.0050600 32 Novel WB029
Lau et al. 2014 WB033 N G 60 4 Arabidopsis thaliana half-sib plant WB033.G WB033.N -0.6590318 0.0000101 33 Novel WB033
Lau et al. 2014 WB033a N G 60 4 Arabidopsis thaliana half-sib plant WB033a.G WB033a.N 1.7678985 0.0000003 34 Novel WB033
Messina et al. 2013 WB035 S G 432 3 Callosobruchus maculatus half-sib insect WB035.G WB035.S -3.0488950 0.2429454 35 Novel WB035
Paccard et al. 2013 WB036 N G 22 7 Arabidopsis lyrata full-sib plant WB036.G WB036.N -0.4328543 0.0080883 36 Novel WB036
Paccard et al. 2013 WB036A N G 22 7 Arabidopsis lyrata full-sib plant WB036A.G WB036A.N -1.1388280 0.0235012 37 Novel WB036
Rauter&Moore 2002 WB038 S G 135 4 Nicrophorus pustulatus half-sib insect WB038.G WB038.S -2.0022426 0.0314474 38 Novel WB038
Relyea 2005 WB039 S G 21 9 Rana sylvatica half-sib amphibian WB039.G WB039.S 0.5733024 0.5624914 39 Novel WB039
Sherrard et al. 2009 WB044 S G 26 10 Avena barbata half-sib plant WB044.G WB044.S 3.4875131 0.5825224 40 Novel WB044
Simons & Roff 1996 WB045 N G 39 7 Gryllus pennsylvanicus full-sib insect WB045.G WB045.N -0.3680302 0.0009672 41 Novel WB045
Simons & Roff 1996 WB045a N G 39 6 Gryllus pennsylvanicus full-sib insect WB045a.G WB045a.N -1.9588013 0.0230411 42 Novel WB045
Simonsen & Stinchcombe 2001 WB046 S G 50 3 Ipomoea hederacea half-sib plant WB046.G WB046.S -2.3591948 0.2837434 43 Novel WB046
Via & Conner 1995 WB053 N G 24 2 Tribolium castaneum half-sib insect WB053.G WB053.N -0.3750967 0.0122579 44 Novel WB053
Via & Conner 1995 WB053_3 N G 24 2 Tribolium castaneum half-sib insect WB053_3.G WB053_3.N 0.3731902 0.1313836 45 Novel WB053
Via & Conner 1995 WB053b N G 80 2 Tribolium castaneum half-sib insect WB053b.G WB053b.N -0.4465280 0.0040389 46 Novel WB053
Via & Conner 1995 WB053b_3 N G 80 2 Tribolium castaneum half-sib insect WB053b_3.G WB053b_3.N 0.9812257 0.0526150 47 Novel WB053
Service & Rose 1985 WB056 N G 64 2 Drosophila melanogaster half-sib insect WB056.G WB056.N -0.7985585 0.1836581 48 Novel WB056
Grill et al. 1997 WB057 N G 26 4 Harmonia axyridis full-sib insect WB057.G WB057.N 1.0197897 0.0237647 49 Novel WB057
Beacham 1988 RR062 A G 10 3 Oncorhynchus gorbuscha half-sib fish RR062.G RR062.A 2.2968023 1.3286533 50 Non-novel RR062
Beacham 1988 RR062_b A G 10 3 Oncorhynchus gorbuscha half-sib fish RR062_b.G RR062_b.A 1.2979916 1.4852452 51 Non-novel RR062
Beacham 1988 RR062b A G 10 3 Oncorhynchus keta half-sib fish RR062b.G RR062b.A 2.4060164 0.8444542 52 Non-novel RR062
Beacham 1988 RR062b_b A G 10 3 Oncorhynchus keta half-sib fish RR062b_b.G RR062b_b.A 0.5874297 0.4801196 53 Non-novel RR062
Evans et al. 2015 RR073 B G 68 8 Poecilia reticulata half-sib fish RR073.G RR073.B -0.4779735 0.2397198 54 Non-novel RR073
Guan et al. 2016 RR076 B G 69 3 Scophthalmus maximus half-sib fish RR076.G RR076.B -0.6184308 0.0032563 55 Non-novel RR076
Sae-Lim et al. 2013 RR087a A G 100 3 Oncorhynchus mykiss full-sib fish RR087a.G RR087a.A -0.4300231 0.0797760 56 Non-novel RR087
Sae-Lim et al. 2013 RR087b A G 100 3 Oncorhynchus mykiss full-sib fish RR087b.G RR087b.A 0.2516590 0.0208090 57 Non-novel RR087
Sae-Lim et al. 2013 RR087c A G 100 3 Oncorhynchus mykiss full-sib fish RR087c.G RR087c.A 0.1919137 0.0845678 58 Non-novel RR087
Kasule 1991 sR097 B G 26 2 Dysdercus fasciatus half-sib insect sR097.G sR097.B -2.4943512 0.4698627 59 Non-novel sR097
Blanckenhorn & Heyland 2004 sR110b B G 37 2 Scathophaga stercoraria half-sib insect sR110b.G sR110b.B 0.0981276 0.0027588 60 Non-novel sR110
Blanckenhorn & Heyland 2004 sR110c B G 40 2 Scathophaga stercoraria half-sib insect sR110c.G sR110c.B -0.3332372 0.0152542 61 Non-novel sR110
Conner et al. 2003 sR115 A G 350 6 Raphanus raphanistrum half-sib plant sR115.G sR115.A -1.8909220 0.0593847 62 Non-novel sR115
Delcourt et al. 2010 sR116 A G 366 8 Drosophila serrata half-sib insect sR116.G sR116.A 3.2396996 0.0071285 63 Non-novel sR116
Tucic et al. 1991 sR145 B G 52 8 Acanthoscelides obtectus full-sib insect sR145.G sR145.B 0.4527517 0.8836976 64 Non-novel sR145
Auld 2010 WB001 B G 30 3 Physa acuta full-sib mollusc WB001.G WB001.B 2.6337736 1.1432721 65 Non-novel WB001
Auld 2010 WB001_b B G 30 3 Physa acuta full-sib mollusc WB001_b.G WB001_b.B -0.9856240 0.2666681 66 Non-novel WB001
Auld 2010 WB001_c B G 30 3 Physa acuta full-sib mollusc WB001_c.G WB001_c.B -0.2218401 0.3625394 67 Non-novel WB001
Begin & Roff 2001 WB002a A G 39 5 Gryllus pennsylvanicus full-sib insect WB002a.G WB002a.A -0.5221598 0.0250811 68 Non-novel WB002
Begin et al. 2004 WB003 A G 56 5 Gryllus firmus full-sib insect WB003.G WB003.A 2.7580171 0.1227235 69 Non-novel WB003
Begin et al. 2004 WB003_2 A G 56 5 Gryllus firmus full-sib insect WB003_2.G WB003_2.A 2.1443823 0.0515381 70 Non-novel WB003
Begin et al. 2004 WB003a A G 45 5 Gryllus firmus full-sib insect WB003a.G WB003a.A -0.1690384 0.0148521 71 Non-novel WB003
Begin et al. 2004 WB003a_2 A G 45 5 Gryllus firmus full-sib insect WB003a_2.G WB003a_2.A 0.4490378 0.0114060 72 Non-novel WB003
Brock et al. 2010 WB007 B G 132 11 Brassica rapa half-sib plant WB007.G WB007.B 0.9015107 0.0265064 73 Non-novel WB007
Collins et al. 1999 WB010 A G 141 4 Achroia grisella half-sib insect WB010.G WB010.A -0.8776703 0.0007712 74 Non-novel WB010
Czesak & Fox 2003 WB011 B G 404 2 Stator limbatus half-sib insect WB011.G WB011.B 1.8636743 0.0002255 75 Non-novel WB011
Delcourt et al. 2009 WB014 A G 319 2 Drosophila serrata half-sib insect WB014.G WB014.A -3.0175489 0.0004793 76 Non-novel WB014
Engqvist 2007 WB016 B G 27 2 Panorpa cognata full-sib insect WB016.G WB016.B -6.7136399 1.2335992 77 Non-novel WB016
Guntrip et al. 1997 WB022 A G 40 2 Callosobruchus maculatus half-sib insect WB022.G WB022.A 1.3420900 0.0578258 78 Non-novel WB022
Holloway et al. 1990 WB025 A G 65 4 Sitophilus oryzae half-sib insect WB025.G WB025.A -1.3078633 0.0037445 79 Non-novel WB025
Klause & Morin 2001 WB028 B G 60 2 Priophorus pallipes full-sib insect WB028.G WB028.B 2.1761862 0.4162618 80 Non-novel WB028
King et al. 2011 WB029 B G 63 3 Gryllus firmus half-sib insect WB029.G WB029.B 1.8019485 0.0051812 81 Non-novel WB029
Lau et al. 2014 WB033 A G 60 4 Arabidopsis thaliana half-sib plant WB033.G WB033.A -0.6597876 0.0000088 82 Non-novel WB033
Lau et al. 2014 WB033a A G 60 4 Arabidopsis thaliana half-sib plant WB033a.G WB033a.A 1.7648037 0.0000004 83 Non-novel WB033
Messina et al. 2013 WB035 B G 432 3 Callosobruchus maculatus half-sib insect WB035.G WB035.B -0.5361305 0.0008421 84 Non-novel WB035
Paccard et al. 2013 WB036 A G 22 7 Arabidopsis lyrata full-sib plant WB036.G WB036.A -0.3174122 0.1071421 85 Non-novel WB036
Paccard et al. 2013 WB036A A G 22 7 Arabidopsis lyrata full-sib plant WB036A.G WB036A.A -1.4644314 0.0091210 86 Non-novel WB036
Rauter&Moore 2002 WB038 B G 135 4 Nicrophorus pustulatus half-sib insect WB038.G WB038.B -2.1617375 0.0489554 87 Non-novel WB038
Relyea 2005 WB039 B G 21 9 Rana sylvatica half-sib amphibian WB039.G WB039.B 0.7823934 0.6700404 88 Non-novel WB039
Sherrard et al. 2009 WB044 B G 26 10 Avena barbata half-sib plant WB044.G WB044.B 4.1801808 1.0238365 89 Non-novel WB044
Simons & Roff 1996 WB045 A G 39 7 Gryllus pennsylvanicus full-sib insect WB045.G WB045.A -0.0302778 0.0128676 90 Non-novel WB045
Simons & Roff 1996 WB045a A G 39 6 Gryllus pennsylvanicus full-sib insect WB045a.G WB045a.A -1.9201388 0.0779041 91 Non-novel WB045
Simonsen & Stinchcombe 2001 WB046 B G 50 3 Ipomoea hederacea half-sib plant WB046.G WB046.B -1.0651373 0.1040128 92 Non-novel WB046
Via & Conner 1995 WB053 A G 24 2 Tribolium castaneum half-sib insect WB053.G WB053.A -0.6445169 0.0009036 93 Non-novel WB053
Via & Conner 1995 WB053_3 A G 24 2 Tribolium castaneum half-sib insect WB053_3.G WB053_3.A 0.7256621 0.0007313 94 Non-novel WB053
Via & Conner 1995 WB053b A G 80 2 Tribolium castaneum half-sib insect WB053b.G WB053b.A -0.6559237 0.0002504 95 Non-novel WB053
Via & Conner 1995 WB053b_3 A G 80 2 Tribolium castaneum half-sib insect WB053b_3.G WB053b_3.A 0.7152326 0.0002055 96 Non-novel WB053
Service & Rose 1985 WB056 A G 64 2 Drosophila melanogaster half-sib insect WB056.G WB056.A -2.0719360 0.0170922 97 Non-novel WB056
Grill et al. 1997 WB057 A G 22 4 Harmonia axyridis full-sib insect WB057.G WB057.A 0.4428487 0.0111591 98 Non-novel WB057

We can now estimate heterogeneity using a multi-level meta-analytic model:

    mod_evolve_angle_int <-  rma.mv(ES ~ 1, V = evolve_angle_data$S_var, random = list(~1|stdy, ~1|err), data = evolve_angle_data)

        stdy_evolve_angle_int <- mod_evolve_angle_int$sigma2[1]
         err_evolve_angle_int <- mod_evolve_angle_int$sigma2[2]

            # Calculate total sampling variance
            wi_evolve_angle <- 1 / evolve_angle_data$S_var
            vi_volve_angle <- sum(wi_evolve_angle) * ( nrow(evolve_angle_data) - 1) / ( (sum(wi_evolve_angle)^2) - sum(wi_evolve_angle))

            I2total_evolve_max_data <- round((stdy_evolve_angle_int + err_sig_angle_within) / (stdy_evolve_angle_int + err_evolve_angle_int + vi_volve_angle), digits = 2)

            I2stdy_evolve_max_data <- round((stdy_evolve_angle_int) / (stdy_evolve_angle_int + err_evolve_angle_int + vi_volve_angle), digits = 2)

But again, we need to control for a number of factors likely to drive variation in effects.

    # Multi-level meta-regression models
        mod_evolve_angle_trait_env <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1 + environment|stdy, ~1|err), data = evolve_angle_data)
        mod_evolve_angle_trait_env_stress <-  rma.mv(ES ~ 1 + environment + design + stress + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1 + environment|stdy, ~1|err), data = evolve_angle_data)
        mod_evolve_angle_trait_env_stress_justNovel <-  rma.mv(ES ~ 1 + design + stress + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1 |stdy, ~1|err), data = subset(evolve_angle_data, environment == "Novel"))

    # calculate marginal mean for non-novel, marginalised over design type. First, refit model
            mod_evolve_angle_trait_env_novel <-  rma.mv(ES ~ 1 + relevel(environment, ref = "Novel") + design + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_angle_data)            
            
            marinal_means_evolve_angle_A_N <- rbind(marginal(mod_evolve_angle_trait_env, "design", data = evolve_angle_data), marginal(mod_evolve_angle_trait_env_novel, "design", data = evolve_angle_data))
            n_A_N <- table(evolve_angle_data$environment)
            marinal_means_evolve_angle_A_N <- cbind(marinal_means_evolve_angle_A_N, n_A_N)
            
        # Now we can do this over design
            mod_evolve_angle_trait_env_half <-  rma.mv(ES ~ 1 + environment + relevel(design, ref = "half-sib") + scale(trait_num), V = evolve_angle_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_angle_data)
            marinal_means_evolve_angle_design <- rbind(marginal(mod_evolve_angle_trait_env, "environment", data = evolve_angle_data), marginal(mod_evolve_angle_trait_env_half, "environment", data = evolve_angle_data))
            n_design <- table(evolve_angle_data$design)
            marinal_means_evolve_angle_design <- cbind(marinal_means_evolve_angle_design, n_design)

            marg_evolve_angle <- rbind(marinal_means_evolve_angle_A_N, marinal_means_evolve_angle_design)

Once we have these models run we can have a look at the table of estimates:

    table_angle_evolve <- data.frame(Est. = round(mod_evolve_angle_trait_env$b, digits =3), L95CI = round(mod_evolve_angle_trait_env$ci.lb, digits=3), U95CI = round(mod_evolve_angle_trait_env$ci.ub, digits =3))
    rownames(table_angle_evolve) <- c("Intercept", "Novel Environment", "Half-sib Breeding Design", "Trait Number")
    kable(table_angle_evolve) %>% kable_styling(font = 12, full_width = TRUE)
Est. L95CI U95CI
Intercept -0.522 -1.566 0.522
Novel Environment -0.105 -0.461 0.251
Half-sib Breeding Design 0.450 -0.788 1.688
Trait Number 0.735 0.215 1.254

5.0 – Supplemental Figures

5.1 – Figure S1 – PRISMA

Figure S1 – PRISMA diagram for the systematic search.

5.2 – Figure S2 – Funnel plots

We explored publication bias using funnel plots, but there was no clear evidence across the effect sizes for publication bias:

Figure S2 – Funnel plots of sampling standard error and effect size. (A-C) ‘Gray’ points represent effect sizes for P, whereas ‘black’ points represent effect sizes for G. (D–F) ‘Orange’ points are effect sizes from novel environments, whereas ‘green’ points are effect sizes from non-novel environments. Note that for C) there was a significant change in angle across environments that was impacted also by the type of environment (stress or not). See main manuscript for more details.

5.3 – Code for Figure 2 - 5 – Main Paper

Figure 2

    #pdf(file = "Fig2.pdf", height = 7.055555, width = 17.513889)
    par(mfrow = c(1,2))

    # Just in case these slots are already in, lets just remove so there are no errors
    varalong_marg_means$name = NULL
    varalong_marg_means$obs = NULL

    marg_list_M1 <- list(SMD_marg_table, varalong_marg_means, angle_marg_means)
    marg_list_M1 <- lapply(marg_list_M1, function(x) cbind(x, name = c("Non-stressful", "Stressful", "Full-sib", "Half-sib")))
    marg_list_M1 <- lapply(marg_list_M1, function(x) cbind(x, obs = c(5,4,1,2)))

    labels = c( "Breeding Design", "Novel Environment")
    for(i in 1){
        par(mar = c(5,13,2,0))
        plot_marginal(marg_list_M1[[i]], "SDV", ylim = c(0,7), xexpan = c(-1.5, 1.5), cex.lab = 2, cex.axis = 1.5)
        #text(marg_list_M1[[i]]$Freq, x = 2, y = marg_list_M1[[i]]$obs, cex = 2)
        mtext(at = 3, side = 2, labels[1], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 6, side = 2, labels[2], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 7, side = 2, "A)", las = 1, font = 2.5, cex = 2.5, adj = 0.9, padj = -0.8)
        text(y = 7, x = -1, expression(bold("Decreased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1,  "Total Variance", las = 1, font = 1, cex = 1.3)
        text(y = 7, x = 1, expression(bold("Increased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1, "Total Variance", las = 1, font = 1, cex = 1.3)
    }


    for(i in 2){
        par(mar = c(5,0,2,13))
        plot_marginal(marg_list_M1[[i]], expression(PV[max]), ylim = c(0,7), xexpan = c(-1.5, 1.5), cex.lab = 2, cex.axis = 1.5, plotname = FALSE)
        text(marg_list_M1[[i]]$Freq, x = 2, y = marg_list_M1[[i]]$obs, cex = 2)
        text("N", x = 2, y = max(marg_list_M1[[i]]$obs)+1, font = 2, cex = 2)
        mtext(at = 7, side = 2, "B)", las = 1, font = 2.5, cex = 2.5, adj = 0.7, padj = -0.8)
        text(y = 7, x = -1, expression(bold("Decreased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1,  expression("variance along"~bold(g[max])), las = 1, font = 3, cex = 1.3)
        text(y = 7, x = 1, expression(bold("Increased")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1, expression("variance along"~bold(g[max])), las = 1, font = 3, cex = 1.3)
    }

    #dev.off()

Figure 3

    #pdf(file = "Fig3.pdf", height = 7.055555, width = 17.513889)
    par(mfrow = c(1,2))
    marg_list_M2 <- list(angle_within_marg_means, marg_evolve_max, marg_evolve_angle)
    marg_list_M2 <- lapply(marg_list_M2, function(x) cbind(x, name = c("Non-novel", "Novel", "Full-sib", "Half-sib")))
    marg_list_M2 <- lapply(marg_list_M2, function(x) cbind(x, obs = c(4,5,1,2)))


    labels = c("Breeding Design", "Environment Type")

    for(i in 2){
        par(mar = c(5,13,2,0))
        plot_marginal(marg_list_M2[[i]], expression(pi[e]), ylim = c(0,7), xexpan = c(-2.5, 1.5), cex.lab = 2, cex.axis = 1.5)
        #text(marg_list_M2[[i]]$Freq, x = 2, y = marg_list_M2[[i]]$obs, cex = 2)
        mtext(at = 3, side = 2, labels[1], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 6, side = 2, labels[2], las = 1, font = 2.5, cex = 1.5, adj = 1)
        mtext(at = 7, side = 2, "A)", las = 1, font = 2.5, cex = 2.5, adj = 0.9, padj = -0.8)
        text(y = 7, x = -1, expression(bold("Low")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1,  "Genetic Variance", las = 1, font = 1, cex = 1.3)
        text(y = 7, x = 1, expression(bold("High")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1, "Genetic Variance", las = 1, font = 1, cex = 1.3)
    }

    for(i in 3){
        par(mar = c(5,0,2,13))
        plot_marginal(marg_list_M2[[i]], expression(theta[e]), ylim = c(0,7), xexpan = c(-1.5, 2.5), cex.lab = 2, cex.axis = 1.5, plotname = FALSE)
        text(marg_list_M2[[i]]$Freq, x = 2, y = marg_list_M2[[i]]$obs, cex = 2)
        text("N", x = 2, y = max(marg_list_M2[[i]]$obs)+1, font = 2, cex = 2)
        mtext(at = 7, side = 2, "B)", las = 1, font = 2.5, cex = 2.5, adj = 0.7, padj = -0.8)
        text(y = 7, x = -1.2, expression(bold("G More Aligned")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = -1.2, "with plasticity vector", las = 1, font = 1, cex = 1.3)
        text(y = 7, x = 1.1, expression(bold("G Less Aligned")), las = 1, font = 3, cex = 1.3)
        text(y = 6.7, x = 1.1, "with plasticity vector", las = 1, font = 1, cex = 1.3)
    }

    #dev.off()

Figure 4 Code

Note that this is not the complete figure. It was modified in Adobe Illustrator.

            Gn_e <- subset(evolve_max_data, environment == "Novel")
            Ga_e <- subset(evolve_max_data, environment == "Non-novel")

            Gn_angle_e <- subset(evolve_angle_data, environment == "Novel")
            Ga_angle_e <- subset(evolve_angle_data, environment == "Non-novel")

        par(mfrow = c(1,2), mar = c(6,6,1,1))
        adj = -0.15
        padj = 0.70
        #emax
        plot(Gn_e$ES ~ Ga_e$ES, ylab = expression(bolditalic(pi[e])~" Novel"), xlab = expression(italic(pi[e])~" Non-Novel"), pch = 16, cex = 1, ylim = c(-5, 16), xlim = c(-5, 16), cex.lab = 1.5, las = 1)
        abline(a = 0, b = 1, lty = 2, lwd = 1.5)
        abline(h = 0, col = "gray", lty = 2, lwd = 1.5)
        abline(v = 0, col = "gray", lty = 2, lwd = 1.5)
        
        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = Ga_e$ES, y1 = (Gn_e$ES+(1.96*sqrt(Gn_e$S_var))), length = 0, angle = 90)
        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = Ga_e$ES, y1 = (Gn_e$ES-(1.96*sqrt(Gn_e$S_var))), length = 0, angle = 90)

        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = (Ga_e$ES+(1.96*sqrt(Ga_e$S_var))), y1 = Gn_e$ES, length = 0, angle = 90)
        arrows(x0 = Ga_e$ES, y0 = Gn_e$ES, x1 = (Ga_e$ES-(1.96*sqrt(Ga_e$S_var))), y1 = Gn_e$ES, length = 0, angle = 90)
        #text(x = Ga_e$ES+0.2, y = Gn_e$ES+0.2, Gn_e$study, cex = 0.2)
        mtext(expression(bolditalic("A)")), adj = adj, padj = padj,cex = 2)


        adj = -0.18
        padj = 0.70
        # e angle
        plot(Gn_angle_e$ES ~ Ga_angle_e$ES, ylab = expression(bolditalic(theta[e])~" Novel"), xlab = expression(italic(theta[e])~" Non-Novel"), pch = 16, cex = 1, ylim = c(-8.5, 6), xlim = c(-8.5, 6), cex.lab = 1.5, las = 1)
        abline(a = 0, b = 1, lty = 2, lwd = 1.5)
        abline(h = 0, col = "gray", lty = 2, lwd = 1.5)
        abline(v = 0, col = "gray", lty = 2, lwd = 1.5)

        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = Ga_angle_e$ES, y1 = (Gn_angle_e$ES+(1.96*sqrt(Gn_angle_e$S_var))), length = 0, angle = 90)
        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = Ga_angle_e$ES, y1 = (Gn_angle_e$ES-(1.96*sqrt(Gn_angle_e$S_var))), length = 0, angle = 90)

        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = (Ga_angle_e$ES+(1.96*sqrt(Ga_angle_e$S_var))), y1 = Gn_angle_e$ES, length = 0, angle = 90)
        arrows(x0 = Ga_angle_e$ES, y0 = Gn_angle_e$ES, x1 = (Ga_angle_e$ES-(1.96*sqrt(Ga_angle_e$S_var))), y1 = Gn_angle_e$ES, length = 0, angle = 90)

        #text(x = Ga_angle_e$ES+0.2, y = Gn_angle_e$ES+0.2, Gn_angle_e$study, cex = 0.2)
        mtext(expression(bolditalic("B)")), adj = adj, padj = padj, cex = 2)

        # Identify values where Gn_e CI's overlap the CI of the Ga_e


    dev.copy2pdf(file = "figure3_standardised10Mar.pdf", useDingbats=FALSE)

Figure 4 – A) The total amount of genetic variation in the direction of the plastic response vector as a proportion of the maximum amount of variation in any direction, π_e. Smaller values on the axes (i.e., negative numbers on axis) correspond to situations where genetic variation along the plasticity is low (‘Low’ label) whereas larger positive values correspond to situations where genetic variation along the plasticity vector is high (‘High’ label) Insets show schematic multivariate ellipses for the Non-novel (‘green’) and novel (‘orange’). Dashed arrows are the plastic response vector (the difference between centroids of multi-variate phenotype space across environments). Arrow length corresponds to ‘smaller’ (negative numbers on axis) or ‘larger’ (positive numbers on axis) amounts of genetic variation along the plasticity vector for the Non-novel (‘green’ text) and novel (‘orange’ text) environments. B) the The angle between the plastic response vector and g_max (θ_e) in novel and Non-novel environments across all studies. Labels correspond to values on the axes where the plasticity vector does not align with the major axis of genetic variation, g_max , (‘Not Aligned’) or where g_max aligns with the plasticity vector (‘Aligned’) in both the Non-novel (‘green’ text) and novel (‘orange’ text) environments. Insets show schematic angles between the plastic response vector (‘dashed arrow’) and g_max in the Non-novel (‘green’) and novel (‘orange’) environments. Axes are on transformed scale to facilitate comparisons with analyses. Large negative values (e.g., -8) indicate that the plastic response vector and g_max is fully aligned (i.e., ~0 degrees), whereas large positive values (e.g., 8) indicate that the plastic response vector and g_maxis not aligned (i.e., ~90). Values of 0 indicate that the plastic response vector and g_max are 45 degrees to one another. For both A) and B), studies that deviate from the expected 1:1 (dashed line) under the assumption of a constant G are colored, and indicate changes in the alignment between G and the plastic response vector across environments. ‘Green’ effects are those where the plasticity vector and G aligns better in the Non-novel environment relative to the novel environment. In contrast, ‘orange’ effects are effects where the plasticity vector and G aligns better with the novel environment relative to the Non-novel environment.

Figure 5 Code

    pdf(height = 9, width = 8)
    titles <- c("A) Kasule 1991", "B) Blanckenhorn & Heyland 2005", "C) Via & Connor 1995","D) Kause & Morin 2001")  

    # locations of the data in the dscrp and mats lists.
    findmeans <- c(6,8,28,20)
    findmats <- rbind(c(22,21,24,23),c(29,30,31,32),c(110,109,112,111),c(78,77,80,79))

    #par(mfrow=c(2,2))
    for each of the for studies:
    for(t in 1:4){
     
      # extract matrices
      mat <- list()
      mat[[1]] <- mats[[findmats[t,1]]]
      mat[[2]] <- mats[[findmats[t,2]]]
      mat[[3]] <- mats[[findmats[t,3]]]
      mat[[4]] <- mats[[findmats[t,4]]]
      
      # extract mean values
      mns <- array(NA,dim=c(2,4))
      mns[,1] <- dscrp[[findmeans[t]]]$mean[1:2]
      mns[,2] <- dscrp[[findmeans[t]]]$mean[1:2]
      mns[,3] <- dscrp[[findmeans[t]]]$mean[3:4]
      mns[,4] <- dscrp[[findmeans[t]]]$mean[3:4]#
      # colors for each matrix
      cols <- c("#348B08A0","#54AB28A0","#DE8601A0","#FEA621A0")#
      # calculate eigenvalues and vectors
      eig <- lapply(mat,eigen)
      
      # calculate coordinates of vectors
      crs <- lapply(eig,function(x)x$vectors*matrix(sqrt(x$values),nrow=2,ncol=2,byrow =TRUE))
      
      # calculate the plotting range
      xlim <- ylim <- list()
      for(i in 1:4){
        x1 <- mns[,c(i,i)]+1.5*crs[[i]]
        x2 <- mns[,c(i,i)]-1.5*crs[[i]]
        xlim[[i]] <- c(x1[1,],x2[1,])
        ylim[[i]] <- c(x1[2,],x2[2,])  
        }
      xlim <- range(unlist(xlim))
      ylim <- range(unlist(ylim))#
      # draw plotting plane 
      plot(NA,xlim=xlim,ylim=ylim,xlab="trait 1",ylab="trait 2", main=titles[t])#
      # add ellipses and their axes for the G and P-matrices
      for(i in 1:4){
        ell <- car::ellipse(mns[,i], as.matrix(mat[[i]]), draw=FALSE, radius=1, center.pch=16, center.cex=1)
        polygon(ell[,1],ell[,2], col=cols[i])
        for(j in 1:2){
          segments(mns[1,i],mns[2,i],mns[1,i]+crs[[i]][1,j],mns[2,i]+crs[[i]][2,j])
          segments(mns[1,i],mns[2,i],mns[1,i]-crs[[i]][1,j],mns[2,i]-crs[[i]][2,j])    
          }
      }
    ## add plasticity vector
    segments(mns[1,1],mns[2,1],mns[1,3],mns[2,3],lwd=3, lty=2) 
    Arrows(mns[1,1],mns[2,1],mns[1,3],mns[2,3],lwd=1, segment ="false" ,arr.type="triangle", arr.width=0.3 , arr.length=0.4, arr.adj = 1)   
    }
    dev.off()

Figure 5 – Examples of studies (with two trait matrices) with high (A, B) or low (C, D) alignment between the direction of plastic responses and standing genetic variation. Arrows are plastic response vectors; large ellipses are P matrices and smaller ellipses are G matrices. Colors correspond to the Non-novel (i.e., ‘green’) and novel (i.e., ‘orange’) environments. A) Kasule (25) used the cotton stainer bug (Dysdercus fasciatus) to estimate quantitative genetic variation in fecundity and adult size under stressful conditions using a half-sib breeding design. B) Blanckenhorn & Heyland (26) used a full-sib breeding design to estimate quantitative genetic variation for hind limb size and development time in the yellow dung fly (Scathophaga stercoraria) under stressful conditions. C) Via & Conner (27) studied the flour beetle (Tribolium castaneum) using a half-sib breeding design to estimate quantitative genetic variation for pupal weight and development time under novel (not stressful) conditions. D) Kause & Morin (28) studied the sawfly (Priophorus pallipes), and estimated quantitative genetic parameters for development time and body size using a full-sib breeding design (broad-sense estimates) under stressful environmental conditions.

6.0 – Additional Supplemental Tables

6.1 – Table S2

Table S2 – Model estimates for (1) standardized differences in total trait genetic variance (SDV) between Non-novel and novel environments, (2) the change in proportion of variance along g_max (PV_max) between Non-novel and novel environments, and (3) the angle between g_max (θ_G) between Non-novel and novel environments. 95% confidence intervals are provided. Intercept values are the average effect size for a study with an average number of traits conducting a full-sib breeding design in a benign environment. All other parameters are contrasts from this mean. Estimates in boldface indicate significant effects (P < 0.05).

Est. (1) L 95 % U 95 % Est. (2) L 95 % U 95 % Est. (3) L 95 % U 95 %
Intercept 0.623 0.086 1.159 0.465 -0.501 1.431 -1.97 -3.015 -0.925
Half-sib Breeding Design -1.092 -1.749 -0.434 -1.475 -2.676 -0.275 -1.245 -2.551 0.062
Stressful Environment -0.012 -0.634 0.610 0.751 -0.387 1.889 1.51 0.269 2.752
Number of Traits -0.090 -0.358 0.179 0.454 -0.044 0.952 0.89 0.353 1.427

6.2 – Table S3

Table S3 – Estimates of the (1) the angle between g_max and p_max within environments, (2) total amount of genetic variance in the direction of plastic response vector (as a proportion of g_max), π_e, and (3) the angle between the plastic response vector and g_max, θ_e, in in both the novel and Non-novel environments. 95% confidence intervals are provided. Intercept values are the average effect size for a study with an average number of traits conducting a full-sib breeding design in an Non-novel environment. All other parameters are contrasts from this mean. Estimates in boldface indicate significant effects (P < 0.05).

Est. (1) L 95 % U 95 % Est. (2) L 95 % U 95 % Est. (3) L 95 % U 95 %
Intercept -2.083 -3.781 -0.385 1.274 -0.353 2.900 -0.525 -1.569 0.520
Novel Environment 0.094 -0.251 0.439 0.088 -0.420 0.597 -0.102 -0.458 0.254
Half-sib Breeding Design 1.080 -0.995 3.155 -0.756 -2.689 1.177 0.449 -0.790 1.689
Number of Traits 0.080 -1.876 2.036 -1.074 -1.875 -0.273 0.733 0.213 1.254

7.0 – Functions Used Throughout

Here we provide the details of the various functions that we use throughout so that readers can have a look at what they are doing in more detail.

## Comparing G and P matrices Functions

# Write files from a list
    writeFiles <- function(listFiles, dir="./data/", add = ".csv"){
        names <- paste0(dir, names(listFiles), add)
        
        for(i in 1:length(names)){
            write.csv(listFiles[[i]], file = names[i])
        }
    }

# Function will go into the data folder and grab all the relevant matrix files ("matrix")or the descriptive data ("D"). It will return a list of folder names with the files within these contents, ignoring irrelevant folders. 

    files <- function(dir = "./data/", file_pattern = "[WRs][BR]", type = c("D", "matrix")){
        type <- match.arg(type)
         tmp <- paste0(dir, list.files(dir)[grep(file_pattern, list.files(dir))])

        if(type == "D"){
            file_names <- lapply(tmp, function(x) paste0(x, "/",list.files(x)[grep("Data", list.files(x))]))
        }

        if(type == "matrix"){
            file_names <- lapply(tmp, function(x) paste0(x, "/",list.files(x)[grep("n=", list.files(x))]))
        }
        names(file_names) <- tmp
        return(file_names)
    }


# Importing matrices
    
    import_mat <- function(list){
            tmp <- as.character(unlist(list))
            mat <- suppressWarnings(lapply(tmp, function(x) read.csv(file = x, header = TRUE)))
            names(mat) <- tmp

            return(mat)
    }


## Function for extracting the sample sizes from the file names. Argument takes the string name for the various matrices and extracts the sample sizes, returning a numeric vector.

    match_n <- function(string){
            m <- regexpr("=[0-9]*(_|.)", string)
            tmp <- regmatches(string, m)
            tmp2 <- gsub("=","", tmp)
            tmp3 <- as.numeric(gsub("_","", tmp2))
        return(tmp3)
    }

## Function extracts the study names of the folders from each of the matrices   
    extract_matrix_studynames <- function(string, pattern = "[WRs][BR][0-9][0-9][0-9].*/"){
            m <- regexpr(pattern, string)
            tmp <- regmatches(string, m)
            tmp2 <- gsub("/","", tmp)
            return(tmp2)
    }

# Extract Environments
    
    extract_matrix_env <- function(string, pattern = "/[ABNSU].", type = c("full", "reduced")){
            m <- regexpr(pattern, string)
            tmp <- regmatches(string, m)
            tmp2.1 <- gsub("/","", tmp)
            tmp2.2 <- gsub("_","", tmp2.1)
            tmp2.3 <- gsub("[0-9]","", tmp2.2)

            if(type == "full"){
                tmp3 <- gsub("_","", tmp2.1)
                return(tmp3)
            }

            if(type == "reduced"){
                return(tmp2.3)  
            }
    }


# Extract whether matrix is a G or P
    extract_matrix_type <- function(string, pattern = "_[GP].csv"){
                m <- regexpr(pattern, string)
                tmp <- regmatches(string, m)
                tmp2 <- gsub("_","", tmp)
                tmp3 <- gsub(".csv","", tmp2)
                return(tmp3)

    }


    ## Simulation function. This will use Monte Carlo simulations to generate various statistics from a given matrix, assuming the traits follow a multivariate normal distribution. A number of different statistics can be chosen, including the total variance in a matrix Vt or "totalVar", which can be used to then look at the difference between two matrices; the proportion of variation along g/p max "varAlong", or it can even extract the full set of eigen values from a matrix ("eigen") or only the G/P max of a matrix ("max"). This works by simply taking a matrix, with it's corresponding sample size and assuming the traits follow a multivariate normal distribution. TThe matrix is used to simulate 1000 matrices according to this MVN distribution, generating the various statistics from each simulated dataset that follows this matrix structure. From these we can get a fairly robust estimate of sampling variance for each matrix and this entire simulated distribution can be used to generate effect sizes.
# x = means
     sim <- function(matrix, x, n, sims = 1000, cor = FALSE, type = c("totalVar", "varAlong", "raw_matrix"), logit = TRUE, matrix_fix = c("shrinkage", "bend"), mu = rep(0, nrow(matrix))){
              
                if(corpcor::is.positive.definite(matrix) == FALSE){

                    if(matrix_fix == "shrinkage"){
                        matrix <- corpcor::cov.shrink(matrix)
                    }
                    
                    if(matrix_fix == "bend"){
                        matrix <- make.positive.definite(matrix)
                    }       
                }

                    eigens <- list()
                    
                    if(type == "totalVar"){
                        for(i in 1:sims){
                          d <- MASS::mvrnorm(n, mu, matrix)
                          if(cor == FALSE){
                          cov <- cov(d)
                          } else {cov <- cor(d)}
                          eigens[[i]] <- sum(eigen(cov)$values)
                        }
                    }

                    if(type == "varAlong"){
                        for(i in 1:sims){
                          d <- MASS::mvrnorm(n, mu, matrix)
                          if(cor == FALSE){
                          cov <- cov(d)
                          } else {cov <- cor(d)}
                          
                          if(logit){
                            p <- (eigen(cov)$values[1]) / sum(eigen(cov)$values)
                            eigens[[i]] <- log(p / (1-p))       
                          } else{
                            eigens[[i]] <- (eigen(cov)$values[1]) / sum(eigen(cov)$values)      
                          }
                        }
                    }

                    if(type == "raw_matrix"){
                        for(i in 1:sims){
                          d <- MASS::mvrnorm(n, mu, matrix)
                          if(cor == FALSE){
                          cov <- cov(d)
                          } else {cov <- cor(d)}
                          eigens[[i]] <- cov
                        }
                    }

                    if(type == "raw_matrix"){
                        return(eigens)
                    }else{
                        return(do.call(rbind,eigens))
                    }
     }


# Function for matrix bending by Reinder

    make.positive.definite <- function(x, tol=1e-6){
      eig <- eigen(x, symmetric=TRUE)
      rtol <- tol * eig$values[1]
      if(min(eig$values) < rtol) {
        vals <- eig$values
        vals[vals < rtol] <- rtol
        srev <- eig$vectors %*% (vals * t(eig$vectors))
        dimnames(srev) <- dimnames(x)
        return(srev)
      } else {
        return(x)
      }
    }

# Angle between two vectors a and b
    angle <- function(a,b){
      
      angle <- acos( sum(a*b) / ( sqrt(sum(a * a)) * sqrt(sum(b * b)) ) ) * (180/pi)

      angtmp <- ifelse(angle > 90, 180-angle, angle)

      # Logit transformation of the angle to normalise the distribution.
      ang <- log( (angtmp/90) / (1 - (angtmp/90) ) )
      
      return(ang)

    }

## Use function below to calculate differences between effects within a study and one environment between G and P. Note that P is subtracted from G. NOTE::: "angle" can ONLY be used with a simulated effect of "max". p = the proportion of vectors, rounded to compare for angles and similarity. x = the columns from the metadata frame. Stats is the simulation object.
    
#type = "SMD"
#p = 0.5
#environment = "across"
#stats = sim_object
#x = metadata
    
    compare_within_env <- function(x, stats, type = c("raw_diff", "SMD", "angle"), environment = c("within", "across")){
            
            cols <- as.vector(x$col)   

        if(type == "raw_diff"){
            col_compute <- data.frame( SMD = (stats[, cols][,2] - stats[, cols][,1]) )
        }

        if(type == "SMD"){
            col_compute <- data.frame( SMD = (stats[, cols][,2] - stats[, cols][,1]) / ( (stats[, cols][,1] + stats[, cols][,2] ) / 2) )
        }


        if(type == "angle"){
                # NOTE::: this can ONLY be used with simulated matrices
                matrix1 <- stats[, cols][,1]
                matrix2 <- stats[, cols][,2]
        
                eig1 <- lapply(matrix1, function(x) eigen(x)$vectors[,1])
                eig2 <- lapply(matrix2, function(x) eigen(x)$vectors[,1])

                # Compute angle
                angle <- mapply(function(a,b) angle(a,b), a = eig1, b = eig2)   
                
                col_compute <- as.data.frame(angle)
            }

        if(environment == "within"){
            colnames(col_compute) <- unique(x$nam.env)
        } else{
            colnames(col_compute) <- unique(x$nam.type)
        }
        
            return(col_compute)
        
    }


# After simulating and creating effects sizes, create the data frame. Change the type of effect size that you wish to generate using the type argument. The choice of this really depends on the specific simulated product that you get. The choices are: c("raw_diff", "SMD", "angle_krzanowski", "angle", "similarity_krzanowski"). raw_diff is simply the subtraction of TOTAL variance of the matrix and the SMD is the standardised mean difference, where the difference in variance is standardised by the combined total variance of each of the matrices. This is important to make sure that G and P are put on comparative scales for between environment comparisons of each matrix.

create_data <- function(sim_object, type = c("SMD", "angle"), environment = c("within", "across"), plot = c("yes", "no")){
    
              nam <- extract_matrix_studynames(colnames(sim_object))
              env <- extract_matrix_env(colnames(sim_object), type = "reduced")
        type_name <- extract_matrix_type(colnames(sim_object))

    # Create a metadata file to summarize and check that everything so far has gone smoothly and as expected.
        
        metadata <- data.frame( file = colnames(sim_object), 
                                 nam = nam, 
                                 env = env, 
                                type = type_name, 
                                   n = n, 
                                 col = 1:length(colnames(sim_object)), 
                                 stringsAsFactors = FALSE )

            metadata$nam.env  <- paste0(metadata$nam, ".", metadata$env)
            metadata$nam.type <- paste0(metadata$nam, ".", metadata$type)

    ## Calculate the difference between G & P within an environment
    
    ## Make sure that every study is arranged by environment and G/P in the same way so that calculations between columns are done correctly. Or, if interested in the GP comparison across environments then we can arrange slightly differently

    # Split metadata on study and environment
        if(environment == "within"){

            split_within <-  split(metadata, metadata$nam.env, drop = TRUE)
            split_within <-  lapply(split_within, function(x) dplyr::arrange(x, type))

        } else {

            split_within <- split(metadata, metadata$nam.type, drop = TRUE)
            split_within <-  lapply(split_within, function(x) dplyr::arrange(x, env))

        }
        

    # Calculate for all environments 
        
        effect_size <- do.call(cbind, lapply(split_within, function(x) compare_within_env(x, sim_object, type = type, environment = environment)))


    # From the simulated distributions, calculate the mean and the sampling variance (i.e., the variance of the distribution)
           
              mean <- sapply(effect_size, mean)
                Vi <- sapply(effect_size, var)
    
    #Create a dataframe
    
    if(environment == "within"){
        dat   <- data.table(study = gsub("[.][SNAB]", "", names(mean)), 
                            env = gsub(".+[.]", "", names(mean)), 
                            ES = mean, 
                            S_var = Vi )
    } else{
        dat   <- data.table(study = gsub("[.][GP]", "", names(mean)), 
                            type = gsub(".+[.]", "", names(mean)), 
                            ES = mean, 
                            S_var = Vi )
    }
    
    # Final data for differences in g/pmax within an environment
        
    if(environment == "within"){
        data <- dplyr::arrange(dat, study, env)
    } else{
        data <- dplyr::arrange(dat, study, type)
    }

    if(plot == "yes"){
            tmp <- effect_size %>% gather()
            fig <- ggplot(tmp, aes(value)) + 
            geom_histogram() + 
            facet_wrap(~key, scales = 'free_x') + theme(axis.text.x=element_text(size=rel(0.5))) + labs(x = "Effect Size")
            return(fig)
        } else{
            return(data)    
        }

}
    


## Evolvability Specific Functions

evolveability_along_plastic_vector <- function(G_sim, cen.diff, compare = "max"){
        
                colnames(G_sim) <- extract_matrix_studynames(colnames(G_sim))
            names_cen.diff_evolve <- names(cen.diff)

            # check that names match
            if(sum(colnames(G_sim)!=names_cen.diff_evolve)>0) {stop("Sorry, you have different matrices and response vectors") }
            
            evol <- c()
            
            for(j in names_cen.diff_evolve){
                        
                        G <- G_sim[,j]
                    r_vec <- cen.diff[[j]]
                
                    tmp <- do.call("rbind", lapply(G, function(x) project.matrix(x, r_vec, compare = compare)))
                    colnames(tmp) <- j
                    
                    evol <- cbind(evol, tmp)
            }

            return(as.data.frame(evol))
}


evolveability_angle_along_plastic_vector <- function(G_sim, cen.diff){
        
                colnames(G_sim) <- extract_matrix_studynames(colnames(G_sim))
            names_cen.diff_evolve <- names(cen.diff)

            # check that names match
            if(sum(colnames(G_sim)!=names_cen.diff_evolve)>0) {stop("Sorry, you have different matrices and response vectors") }
            
            evol <- c()
            
            for(j in names_cen.diff_evolve){
                
                    # Grab the same study for response vector and G     
                        G <- G_sim[,j]
                    r_vec <- cen.diff[[j]]
                    
                    # Calculate the angle between gmax and response vector
                    tmp <- do.call("rbind", lapply(G, function(x) angle(eigen(x)$vectors[,1], r_vec)))
                    colnames(tmp) <- j
                    
                    evol <- cbind(evol, tmp)
            }

            return(as.data.frame(evol))
}

create_data_evolve <- function(sim_object, meta_evolve, env = c("N", "S")){
    
        columns <- c("authors", "year", "study", "env", "type", "n", "trait_num", "genus", "species", "design", "group", "name.type", "name.env")
    trim_meta <- meta_evolve[meta_evolve$type == "G", columns] 
    trim_meta <- trim_meta[trim_meta$env %in% env,]

            mean <- sapply(sim_object, mean)
             Vi <- sapply(sim_object, var)
             study <- names(mean)

             dat   <- data.table(study=study,  
                            ES = mean, 
                            S_var = Vi )

             data_fin<- left_join(trim_meta, dat, by = "study")

             return(data_fin)
}


combine_meta_data_evolve <- function(es_data, meta_sub){
            
                es_data <- unique(left_join(es_data, meta_sub, by = c("file")))

            # Fix up a couple columns
            es_data$species_full <- paste0(es_data$genus, "_", es_data$species)
            es_data$stdy <- substr(es_data$study, 1, 5)
            rownames(es_data) <- 1:nrow(es_data)
        return(es_data)
}


sim_G <- function(Gmatrix, n, sims = 1000, cor = FALSE, matrix_fix = c("shrinkage", "bend")){
                    
                    mu <- rep(0, nrow(Gmatrix))
                
                if(corpcor::is.positive.definite(Gmatrix) == FALSE){

                    if(matrix_fix == "shrinkage"){
                        Gmatrix <- corpcor::cov.shrink(Gmatrix)
                    }
                    
                    if(matrix_fix == "bend"){
                        Gmatrix <- make.positive.definite(Gmatrix)
                    }       
                }

                matrices <- list()
                    
                        for(i in 1:sims){
                          
                          d <- MASS::mvrnorm(n, mu, Gmatrix)
                          
                          if(cor == FALSE){
                            cov <- cov(d)
                          } else {
                            cov <- cor(d)
                          }
                          
                          matrices[[i]] <- cov
                        }       

                return(matrices)                    
}

centroid.diff <- function(a,b){
  # a: vector a
  # b: vector b
  dz <- b-a
  val <- sqrt(sum(dz^2))
  vec <- dz/val
  return(list(value=val,vector=vec))
}

# Function to project a covariance matrix onto a vector.

project.matrix <- function(mat,vec,compare=NA){
  # mat: covariance matrix that needs to be projected
  # vec: vector on which the covariance matrix needs to be projected
  # compare: to what does the amount of variance of mat on vec needs to be compared to:
  #   NA        : no comparison, output is the amount of variance on vec
  #   "max"     : as proportion of the first eigenvector of mat
  #   "maxcorr" : as proportion of the first eigenvector of mat, corrected for the expected value
  #   "mean"    : as proportional difference to mean value of all eigenvectors of mat
  #   "total"   : as proportion of the sum of the values of all eigenvectors of mat
  
  dims       <- length(vec)
  beta_prime <- matrix(vec,ncol=dims,nrow=dims)
  beta       <- matrix(vec,ncol=dims,nrow=dims,byrow=TRUE)
  abs_beta_2 <- sum(vec^2)
  pv <- as.numeric((t(vec)%*%mat%*%vec)/abs_beta_2)
  ev <- eigen(mat)$values
  if(is.na(compare)){
    return(pv)
  }else if(compare=="max"){
    return(log( (pv/ev[1]) / (1 - (pv/ev[1]))))
  }else if(compare=="maxcorr"){
    return(log( (pv/ev[1]) / (1 - (pv/ev[1])))-log( (mean(ev)/ev[1]) / (1 - (mean(ev)/ev[1]))))
  }else if(compare=="total"){
    return(pv/sum(ev))    
  }else{
    return((pv-mean(ev))/mean(ev))    
  }
}


## Adding some new functions for calculating marginalised means. 

# Relevels the data based on a factor in the model and a reference level of the factor
    re_level_dat <- function(data, factor, ref){
        data[, factor] <-   relevel(data[, factor], ref = ref)
        return(data)
    }

    #@example test <- re_level_dat(evolve_max_data, "design", ref = "half-sib")


## Extract coefficients from models
    extract_coefs <- function(model){
        mean <- model$b[1]
          se <- model$se[1]

          dat <- data.frame(mean = mean, se = se)

     return(dat)
    }

# Calculate marginal means for a single factor
    marg_mean <- function(coefs){
        mean_mar <- ( (coefs[1,"mean"]*coefs[1,"n"]) + (coefs[2,"mean"]*coefs[2,"n"]) ) / sum(coefs$n)
        return(mean_mar)
    }

    #@examples marg_mean(coefs)

#Calculate marginal SE's for a single factor

    marg_se <- function(coefs){
        n1 <- (coefs[1,"n"] / sum(coefs$n))
        n2 <- (coefs[2,"n"] / sum(coefs$n))

        mar_se  <- sqrt((coefs[1,"se"]^2)*(n1^2) + (coefs[2,"se"]^2)*(n2^2))
        return(mar_se)
    }

    #@example marg_se(coefs)

# Function for producing marginalised mean and se across a single factor level. To get estimates for a 2x2 factor design simply refit the model re-leveling the second factor (i.e., not the one used in the "marginal" function)
    
    marginal <- function(model, factor, data){
            
        # Get relevant data
                 ref <- levels(data[ , factor])
            n_levels <- table(data[ , factor])

        # Get first set of estimates from model
               coef1 <- cbind(extract_coefs(model), name = ref[1])

        # Update data, to get the new estimates
             dat_new <- re_level_dat(data, factor, ref = ref[2])

        # Refit model
            model_new <- update(model, data = dat_new)

        # Extract updated coefs
            coef2 <- cbind(extract_coefs(model_new), name = ref[2])

        # Concatenate the two rows
            coefs <- cbind(rbind(coef1, coef2), n = as.vector(n_levels))

        #Calculate the marginal mean and SE across the levels of "factor"
            marginal <- data.frame(mean = marg_mean(coefs),
                                   se   = marg_se(coefs)  )
        
        return(marginal)
    }

 #  @example{
#           full_A <-  rma.mv(ES ~ 1 + environment + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)#

#           marginal(full_A, "design", evolve_max_data) # Marginal mean of Non-novel across design#

#           full_N <-  rma.mv(ES ~ 1 + relevel(environment, ref = "Novel") + design + scale(trait_num), V = evolve_max_data$S_var, random = list(~1+environment|stdy, ~1|err), data = evolve_max_data)#

#           marginal(full_N, "design", evolve_max_data) # Marginal mean of Novel across design
#   }


    plot_marginal <- function(table, effectsize, xexpan, plotname = TRUE, ...){
        plot(obs ~ mean, type = "n", las = 1, ylab = "", xlab = effectsize, data = table, bty = "n", yaxt = "n",  xlim = range(table$mean)+xexpan, ...)
        abline(v=0, lty = 2, lwd = 1.5)
        points(obs ~ mean, data = table, pch = 16, cex = 1.8)

        arrows(x0=table$mean, y0=table$obs, x1=table$mean + 1.96*(table$se), y1=table$obs, length = 0, angle = 90)
        arrows(x0=table$mean, y0=table$obs, x1=table$mean - 1.96*(table$se), y1=table$obs, length = 0, angle = 90)  

        if(plotname){
            mtext(table$name, at = table$obs, side = 2, cex = 1.3, las = 1)
        }
    }